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ABSTRACT 



The Laser Interferometer Gravitational Wave Observatory (LIGO) is one of a 
new generation of detectors of gravitational radiation. The existence of gravitational 
radiation was first predicted by Einstein in 1916, however gravitational waves have 
not yet been directly observed. 

One source of gravitation radiation is binary inspiral. Two compact bodies orbit- 
ing each other, such as a pair of black holes, lose energy to gravitational radiation. 
As the system loses energy the bodies spiral towards each other. This causes their 
orbital speed and the amount of gravitational radiation to increase, producing a char- 
acteristic "chirp" waveform in the LIGO sensitive band. 



iii 



In this thesis, matched filtering of LIGO science data is used to search for low mass 
binary systems in the halo of dark matter surrounding the Milky Way. Observations 
of gravitational microlensing events of stars in the Large Magellanic Cloud suggest 
that some fraction of the dark matter in the halo may be in the form of Massive 
Astrophysical Compact Halo Objects (MACHOs). It has been proposed that low 
mass black holes formed in the early universe may be a component of the MACHO 
population; some fraction of these black hole MACHOs will be in binary systems and 
detectable by LICO. 

The inspiral from a MACHO binary composed of two 0.5 solar mass black holes 
enters the LICO sensitive band around 40 Hz. The chirp signal increases in amplitude 
and frequency, sweeping through the sensitive band to 4400 Hz in 140 seconds. By 
using evidence from microlensing events and theoretical predictions of the population 
an upper limit is placed on the rate of black hole MACHO inspirals in the galactic 
halo. 
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Conventions 



There are two possible sign conventions for the Fourier transform of a time domain 
quantity v{t). In this thesis, we define the Fourier transform v{f) of a v{t) to be 



This convention differs from that used in some gravitational wave literature, but is 
the adopted convention in the LIGO Scientific Collaboration. 

The time-stamps of interferometer data arc measured in Global Positioning System 
(GPS) seconds: seconds since 00:00.00 UTC January 6, 1980 as measured by an 
atomic clock. 

Astronomical distances arc quoted in parsecs 




and the inverse Fourier transform to be 




1 pc = 3.0856775807 x 10 



m 



and masses in units of solar mass 



IM, 



■0 — 



1.98892 X 10' 
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Chapter 1 



Introduction 



One of the earliest predictions of the Theory of General Relativity was the existence 
of gravitational waves. By writing the metric g^i, as the sum of the flat Minkowski 
metric ?7^j^ and a small perturbation /i^j,, 



and considering bodies with negligible self-gravity, Einstein showed [T] 

"that these /i^j^ can be calculated in a manner analogous to that of the 
retarded potentials of electrodynamics." 

It follows that gravitational fields propagate at the speed of light. In electrodynamics, 
the lowest multipole moment that produces radiation is the electric dipole; there is 
no electric monopole radiation due to the conservation of electric charge. Similarly 
in General Relativity, the lowest multipole that produces gravitational waves is the 
quadrupole moment. Radiation from the mass monopole, mass dipole and momen- 
tum dipole vanish due to conservation of mass, momentum and angular momentum 
respectively. Einstein also derived the quadrupole formula for the gravitational wave 
field, which states that the spacetime perturbation is proportional to the second time 
derivative of the quadrupole moment of the source. The strength of the gravitational 
waves decreases as the inverse of the distance to the source. We can estimate this 
strength at a distance r by noticing that the quadrupole moment involves terms of 
dimension mass x length^ and so the second time derivative of the quadrupole mo- 
ment is proportional to the kinetic energy of the source associated with non-spherical 
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motion E^^. Using the quadrupole formula, which we will see in equation (2.63), we 
then approximate the strength of the gravitational waves as 

— ■ (1-2) 

The effect of a gravitational wave is to cause the measured distance L between two 
freely falling bodies to change by a distance AL ~ hL. 

Interferometers were suggested as a way of measuring the change in length between 
two test masses by Pirani in 1956 [2] and the first working detector was constructed 
by Forward in 1971 [3J. The fundamental designs of modern laser interferometers were 
developed by Weissjl] and DreverjH] in the 1970s. The principle upon which interfero- 
metric detectors operate is to use laser light to measure the change in distance between 
two mirrors as a gravitational wave passes through the detector. The sensitivity of an 
interferometer on the Earth is limited by gravity gradient noise at frequencies below 
~ 5 Hz[n]. Any time changing distribution of matter near the detector, for example 
compression waves in the Earth, cause fluctuations in the local gravitational field. 
These fluctuations will cause the test masses to move producing a spurious response 
in the interferometer which masks the presence of gravitational waves. In fact. Earth 
based interferometers are typically limited in sensitivity to frequencies above ~ 10 Hz 
due to the seismic motion of the earth. 

The canonical example of an astrophysical source of gravitational waves is the 
Hulse- Taylor binary pulsar, PSR 1913+ 16 [7]. This system is composed of two neutron 
stars, each of mass ~ 1.4 Mq, with average separation and orbital velocity of ~ 10^ m 
and ~ lO^ms"^, respectively. The period of the orbit is 7.75 hours and the binary is 
at a distance from the earth of ~ 6 kpc. Hulse and Taylor observed that the orbital 
period of the binary is decreasing and that the rate of orbital energy loss agrees 
with the expected loss of energy due to the radiation of gravitational waves to within 
0.3%[Sl IS]- Since the quadrupole moment of an equal mass binary is periodic at half 
the orbital period, we would expect the frequency of the gravitational waves emitted to 
be twice the orbital frequency. Thus, the gravitational waves from PSR 1913 + 16 have 
a frequency /gw ~ 10""^ Hz that is outside the sensitive band of earth based detectors. 
Nevertheless, the orbit will continue to tighten by gravitational wave emission, and 
the two neutron stars are expected to merge in about 300 million years; in the last 
several minutes prior to merger, the gravitational wave frequency will sweep upward 
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from ~ 10 Hz reaching about 1500 Hz just before the merger. 

It is worthwhile to estimate the strength of the gravitational waves from a neutron 
star binary since it informs the target sensitivity for modern interferometric detectors. 
The non-spherical kinetic energy of this system is 

KL~1.4MQ(7ra/Tf (1.3) 

where T is the binary period and a is the average separation. The period, separation 
and mass of a binary are related by Kepler's third law, 



where M is the total mass of the binary. Using equation (1.2), Kepler's third law and 



the non-spherical kinetic energy given in equation (1.3), we can estimate the strength 
of the waves from a neutron star binary as 



When the orbital separation is a ~ 10^ m, the orbital period will be T ~ 10~^ seconds 
and the gravitational wave strain will be /i ~ 10^^^. 

To date, four binary neutron star systems that will merge within a Hubble time 
have been discovered. By considering the time to merger, position and efficiency of 
detecting such binary pulsar systems, the galactic merger rate for inspirals can be 
estimated [lOj. The latest estimates of neutron star inspirals in the Milky Way are 
8.3 X 10~^ yr~^- Extrapolating this rate to the neighboring Universe using the blue- 
light luminosity gives an (optimistic) estimate of the rate at 0.3 yr~^ within a distance 
of ~ 20 Mpc. To measure the waves from a neutron star binary at this distance, we 
must construct interferometers that are sensitive to gravitational waves of strength 
h ~ 10~^^. An overview of the theory and experimental techniques underlying the 
generation and detection of gravitational waves from binary inspiral is presented in 
chapter |2j 

A world-wide network of gravitational wave interferometers has been constructed 
that have the sensitivity necessary to detect the gravitational waves from astrophysical 
sources. Among these is the Laser Interferometric Gravitational Wave Observatory 
(LIGO)[TT]. LIGO has completed three science data taking runs. The first, referred 



to as SI, lasted for 17 days between August 23 and September 9, 2002; the second, 
S2, lasted for 59 days between February 14 and April 14, 2003; the third, S3, lasted 
for 70 days between October 31, 2003 and January 9, 2004. During the runs, all 
three LIGO detectors were operated: two detectors at the LIGO Hanford observatory 
(LHO) and one at the LIGO Livingston observatory (LLO). The detectors are not 
yet at their design sensitivity, but the detector sensitivity and amount of usable data 
has improved between each data taking run. The noise level is low enough that 
searches for coalescing compact neutron stars are worthwhile, and since the start of 
S2, these searches are sensitive to extra-galactic sources. Using the techniques of 
matched filtering described in chapter |4] of this dissertation, the SI binary neutron 
star search set an upper limit of 

7^9o% < 1-7 X 10^ per year per Milky Way Equivalent Galaxy (MWEG) (1.6) 

with no gravitational wave signals detected. Details of this analysis can be found in 

In this dissertation, we are concerned with the search for gravitational waves from 
a different class of compact binary inspiral: those from binary black holes in the 
galactic halo. Observations of the gravitational microlensing of stars in the Large 
Magellanic cloud suggest that ~ 20% of the galactic halo consists of objects of mass 
~ 0.5 of unknown origin. In chapter [3] we discuss a proposal that these Massive 
Astrophysical Compact Halo Objects (MACHOs) may be black holes formed in the 
early universe and that some fraction of them may be in binaries whose inspiral is 
detectable by LIGO [13]. The upper bound on the rate of such binary black hole 
MACHO inspirals are projected to be i? ~ 0.1 yr^^ for initial LIGO, much higher 
than the binary neutron star rates discussed above. It should be noted however, 
that while binary neutron stars have been observed, there is no direct observational 
evidence of the existence of binary black hole MACHOs. Despite this, the large 
projected rates make them a tempting source for LIGO. In chapter |5] we describe an 
analysis pipeline that has been used to search the LIGO S2 data for binary black hole 
MACHO^ Chapter [6] describes how the search techniques were tested on data from 
the gravitational wave interferometers. Finally we present the result of the S2 binary 

"'^The same pipeline has also been used to search for binary neutron star inspiral in the S2 data 
and the results of this search will be presented in jT4] . 
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black hole MACHO search in chapter [7} 
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Chapter 2 

Gravitational Radiation from 
Binary Inspiral 

In this chapter we review some of tlie pliysics underlying the detection of gravitational 



waves from binary inspiral. In section 2.1| we review the effect of gravitational waves 



on a pair of freely falling particles in order to introduce some of the concepts that 
we need to discuss the detection of gravitational waves from binary inspiral. For a 



detailed description of gravitational waves, we refer to [151 [IS]- Section 2.2 describes 
how a laser interferometer can be used to measure this effect. The gravitational 
waveform produced by the inspiral of two compact objects, such as neutron stars or 



black holes, are discussed in section |2.3 We also derive the waveform that will be 



used to search for gravitational waves from binary inspiral events in the Universe. 

2.1 The Effect of Gravitational Waves on Freely Falling Par- 
ticles 

The 4- velocity m of a freely falling test particle satisfies the geodesic equation [T7j 

(Vs^r)° = = 0, (2.1) 

where ; denotes the covariant derivative, that is, 

u-.X = + V) (2-2) 

where F"^ is the connection coefficient of the metric Qf^^ and , fi represents the standard 
partial derivative with respect to the coordinate x^. 
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Consider two particles A and B, as shown in figure [T] (a), with separation vector 
^. The particles are initially at rest with respect to each other, so 

V^u = 0, (2.3) 
M-e = 0. (2.4) 

If the spacetime is curved, the second derivative of ^ along u is non-zero; it is given 
by the equation of geodesic deviation 

V„^V„^f=-i?(_,u,e,M), (2.5) 

where i?(_, u, C,, u) is the Riemann curvature tensor. If the spacetime is flat with weak 
gravitational waves propagating in it, we can describe it by a metric 

g^y = Vf^u + h^u, (2.6) 

where h^j^^, is the perturbation to the metric due to the gravitational waves and r^^j^ 
is the flat Minkowski metric . We now introduce a Local Lorentz Frame (LLF) for 
particle A. The LLF of particle A is a coordinate system x° in which 

9^.u{A) = ri^, (2.7) 

and 

g^^uA^) = 0, (2.8) 

where g^i^{A) is the value of the metric at point A. This LLF is equivalent to a 
Cartesian coordinate system defined by three orthogonally pointing gyroscopes carried 
by particle A. The curvature of spacetime means that the coordinate system is not 
exactly Cartesian, but it can be shown that this deviation is second order in the 
spatial distance from the particlep^. This means that along the worldline of particle 
A the metric is 

g,u = ii^u + (2.9) 

where x is the distance from the particle and R ~ \Rai3'y5\- We can write the Cartesian 
coordinates of the LLF of A a,s = {x^ , x*), where x° = t is the timelike coordinate 
and X* are the three Cartesian coordinates. Then in the LLF of particle A the equation 
of geodesic deviation becomes 

^ = -R\p,n'^eu' = -R\,o^\ (2.10) 



8 



since u = (1,0,0,0). The presence of the gravitational waves are encoded in the 
curvature Rap-yS which satisfies the wave equation 

V^''RaP^5,^.u = 0. (2.11) 

In the Local Lorentz frame, the components of ^ are just the coordinates of B. In 
the LLF of A we may write 

e = e?o) + '^e\ (2.12) 

where ^^^^ is the unperturbed location of particle B and 5^^ is the change in the 



position of B caused by the gravitational wave. Substituting equation (2.12) into the 
equation of geodesic deviation, we obtain 



d'^s^j fc k 

Q^2 ~ ~^ 0A:0^(0) = --RjOfcO^(0)5 (2.13) 



where we have used r/^i, to lower the spatial index j of the Riemann tensor. For a weak 
gravitational wave, all the components of Rap-yS are completely determined by Rjow- 
Furthermore, it can be shown that the 3x3 symmetric matrix Rjoko, which we would 
expect to have 6 independent components, has only 2 independent components due to 
the Einstein equations and the Biancci identity. We define the (transverse traceless) 
gravitational wave field, hj^, by 



Using this definition in equation (2.13), we obtain 

= Ihji^^f.y (2.15) 

If we orient our coordinates so the gravitational waves propagate in the ^-direction, 

so hjl^{t — z), then the only non-zero components of hj^ are /iJJ, h^y, h^J and hyj. 
Since hj^ is symmetric and traceless, these components satisfy 

hll = -h^y, (2.16) 

h!y = h^I- (2.17) 

For two more particles C and D separated by 



(2.18) 
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as shown in figure [T] (b), the effect of the gravitational wave is then given by 

= IhJ^Ctoy (2.19) 

Taking the two particles A and B to he on the x-axis of the LLF of particle A with 
separation X(o), without loss of generality, we may write 

^=(x(o) + 5x,0,0), (2.20) 

where 5x is the displacement of particle B caused by the gravitational wave. Similarly, 
if particles C and D lie on the y-axis of the LLF of particle C with separation ?/(o), 
we may write 

C = (0,i/(o) + <5y,0) (2.21) 

where 6y is the displacement of particle D caused by the gravitational wave. 
We define the two independent components of the gravitational wave to be 

h+ = h^I = -h^yl (2.22) 
K = h^J = h^J (2.23) 

which we call the plus and cross polarizations of the gravitational wave respectively. 
The influence of a linearly + polarized gravitational wave propagating in the z- 



direction on the particles A,B,C,D is then given by substituting equations (2.20) 



and (2.21) into (2.15) and (2.19) respectively to obtain 

1 



Sx{t-z) = -h^J{t-z)x^o), (2.24) 

5yit-z) = -h^;^it-z)y^o) (2.25) 

Similarly, for a linearly x polarized gravitational wave propagating in the z-direction 
the effect on the particles is 

1 

h 

2 ^'^ 

6y{t-z) = lh^J{t-z)x(^o). (2.27) 



Sx{t-z) = ^h]:^{t-z)y^o), (2.26) 
2"yx \" ~/"^(0). 



Figure |2] shows the effect of and h^. on a. ring of particles that lie in the xy plane. 
We can see for the plus polarization that the effect of a gravitational wave is to stretch 
the ring in the x direction, while squeezing it in the y direction for the first half of a 
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cycle and then squeeze in the x direction and stretch in the y direction for latter half 
of the cycle. There is therefore a relative change in length between the two particles 
AB and CD as a gravitational wave passes. The overall effect of a gravitational wave 
containing both polarizations propagating in the z direction is 

6x{t ~z) = \ {K{i - ^)a;(o) + /^x (t - 2:)y(o)] , (2.28) 
by(t -A = \ {-K[t - z)y^o) + K{t- 2)x(o)] . (2.29) 
It is the change in the distance between a pair of particles that we attempt to measure 



with gravitational wave detectors. We can see from equation (2.29) that the change 
in length is proportional to the original distance between the test masses. For a pair 
of test masses separated by a length L, we define the gravitational wave strain h to 
be the fractional change in length between the masses 

h^--. (2.30) 

The reason to include a factor of 1 /2 in this definition will become apparent when we 
discuss measuring gravitational wave strain with an interferometer. 



2.2 The LIGO Gravitational Wave Detectors 

Several major efforts are underway [TT| [T8l 119] to measure the strain produced by a 
gravitational wave using laser interferometry. The results in this thesis are based on 
data from the Laser Interferometer Gravitational wave Observatory (LIGO). LIGO 
operates three power-recycled- Fabry- Perot-Michels on interferometers in the United 
States. Two of these are co-located at the LIGO Hanford Observatory, WA (LHO) 
and one at the LIGO Livingston Observatory (LLO). The interferometers at LHO 
are 4 km and 2 km in arm length and are referred to as HI and H2, respectively. 
The interferometer at LLO is a 4 km long interferometer referred to as LI. The 
locations and names of the detectors are shown in figure |3} As we saw in chapter [T| 
to detect the gravitational wave strain produced by typical astrophysical sources we 
need to measure h ~ 10~^^. If we separate our test masses by a distance of 4 km (a 
practical distance for earthbound observatories) the challenge faced by gravitational 
wave astronomers is to measure changes of length of order 

AL ~ 10^2^ X lO^m ~ lO^^^m. (2.31) 
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2.2.1 The Design of the LIGO Interferometers 

In an interferometric gravitational wave detector the freely falling masses described 
in the previous section are the mirrors that form the arms of the interferometei{^ 
and laser light is used to measure the change in length between the mirrors. The 
challenge facing experimenters constructing a gravitational wave interferometer is to 
measure changes of length of order ~ 10~^^m using laser light has a wavelength of 
A; ~ 10"^ m. It should be noted that measuring a phase shift 

A$ ~ ^ ~ 10-12 (2.32) 

is a factor 10^^ more sensitive than the interferometers used by Michelson and Morely 
to disprove the existence of the ether. 

A schematic of a simple Michelson interferometer is illustrated in figure |4] (a). 
Laser light is shone on a beam splitter which refiects half the light into the x-arm 
and transmits half the light into the y-arm of the interferometer. The light travels a 
distance L in each arm and then is refiected back towards the beam splitter by the end 



test masses. These masses are equivalent to the test masses B and D in section |2.1 
Consider the light in the x-arm. For the laser light, the spacetime interval between 
the beam splitter and the end test mass is given by 

ds"^ = g^^dx'^dx'' = 0. (2.33) 

In the presence of a plus polarized, sinusoidal, gravitational wave traveling in the 



2;-direction, equation (2.33) becomes 



c'^dt'^ = [l + h+{t- z)] dx"^ + [l-h+{t- z)] dy^ + dz"^. (2.34) 

We can measure the response of the interferometer to a gravitational wave by consider- 
ing the phase shift of light in the arms. The phase that the light acquires propagating 



^Thc mirrors in an Earth bound gravitational wave observatory are not truly freely falling as 
they are accelerated by the gravitational field of the Earth. It can be shown that the horizontal 
motion of suspended mirrors is the same as that of freely falling test masses. 
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from the beam splitter to the x-end test mass and back is given by[6] 



= / 271 f I dt 







If If 

-/ 2Tifi^l + h+dx- - / 2tx fi^/l^hl dx (2.35) 



c 

where trt is the round trip time of the hght and // is its frequency. We have discarded 
higher order terms in hj^ as their effect is neghgible. We can see that the phase shift 
acquired in the x-arm due to the gravitational wave is 

5$^ = '^h+L. (2.36) 

A similar calculation shows that the phase shift acquired in the y-arm is 

27r 

5^y = - — h+L (2.37) 

and so the difference in phase shift between the arms is 

A<l> = ^h+L. (2.38) 

A typical astrophysical source of gravitational radiation of interest to LIGO, has a 
frequency /gw ~ 100 Hz. Therefore the wavelength of the gravitational wave is 
Aqw ~ 3000 km. If trt = 1//gw there will be no phase shift of the light at leading 
order in h^. The light spends exactly one gravitational wave period in the arm and so 
the phase shift acquired by positive values of h^(t — z) is canceled out by the phase 
shift due to negative values of hj^{t — z). The interferometer achieves maximum 
sensitivity when the light spends half a gravitational wave period in the arms, that is 

L = ^~ 1000 km (2.39) 

which is a hopelessly impractical length for a earthbound detector. Instead, a simple 
Michelson interferometer is enhanced by placing two additional mirrors in the arms 
of the interferometer near the beam splitter, as shown in figure |4] (b). These inner x 
and y test masses (referred to as ITMX and ITMY) are designed in LIGO to store 
the light in the arms for approximately one half of a gravitational wave period. The 
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mirrors create a Fabry-Perot cavity in each arm that stores the hght for B ~ 200 
bounces, giving a phase shift of 

L 4x1 0*^ III 
A$ = 4:7i—Bh+ ~ 10 X X 200 x h+. (2.40) 

\i + 10-6 m + ^ ^ 

For a gravitational wave strain of /i ~ 10^^^, this increases A$ by 3 orders of mag- 
nitude to a phase shift A$ ~ 10~^. Further increasing B does not gain additional 
sensitivity, however, as storing the light for longer than half a gravitational wave pe- 
riod causes it to lose phase shift as the sign of the gravitational wave strain changes. 

Is it possible to measure a phase shift of 10~^ using a Fabry-Perot-Michelson 
interferometer? We measure the phase shift by averaging the light at the photodiode 
over some period, r. Let N be the number of photons from the laser arriving at the 
photodiode in the time r. The measured number of photons in the averaging interval 
is a Poisson process, with probability distribution function for N given by 

PiN) = ^, (2.41) 

where N is the mean number of photons per interval r. The la uncertainty in the 

number of photons arriving in the averaging time is therefore 



AN = VN. (2.42) 

The accuracy to which we can measure the phase shift for a given input laser power 
is constrained by the uncertainty principle, 

AtAE>^ (2.43) 

as follows. The energy of the light arriving at the photodiode in time r is 

E = h^N, (2.44) 

which, due to the counting of photons, has uncertainty 

AE = n^VN. (2.45) 

The uncertainty in the measured the phase is related to the uncertainty in the time 
that a wavefront reaches the beam splitter , i.e. 

At 

A$ = 27rc— (2.46) 

A; 



Substituting equation (2.45) and (2.46) into equation (2.43), we obtain 

AtA£;=^fiH^£^>5 (2.47) 

The accuracy which which we can measure the phase is therefore no better than 

A$ > (2.48) 

Hence photon counting statistics hmits the accuracy with which the phase shift can 
be measured by this method, and this equation tells us how many photons we need 
in an averaging period to measure a given phase shift. We need at least 

> ^ (2.49) 

~2(A<I>)' ^ ^ 

photons to measure the phase shift. The optimal averaging time for a gravitational 
wave with frequency /gw is half a period so that the light acquires the maximum 
phase shift, that is 

TR.-^. (2.50) 

The intensity of laser light required to measure a phase shift of 10~^ for a gravitational 
wave of /gw ~ 100 Hz is then 



27:hc\ / 1 ^ ^ 



■^i J \2/gw, 

2/gw (2.51) 



1 V f2Tlhc 



1 \ VlO^ X 10"34 X 10* 



10^ ~ 10^ w, 



10- V V 10-^ 

however, lasers used in the first generation of interferometers have a typical output 
power of ~ 5 W. To increase the power in the interferometer, the final enhancement 
to the basic design of our interferometer is the addition of a power recycling mirror 
(RM) between the beam splitter and the laser, as shown in figure |4] (c). This mirror 
reflects some of the (otherwise wasted) laser light back into the interferometer and 
increases the power incident on the beam splitter so that the phase shift due to a 
gravitational wave of order h ~ 10^^^ can be measured. 

The laser light must be resonant in the power recycling and Fabry-Perot cavities, 
to achieve the required power build up in the interferometer. This requires a compli- 
cated length sensing and control system^^ which continuously monitors the positions 
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of the mirrors in the interferometer and apphes feedback motions via electromagnetic 
actuators. The interferometer is said to be locked when the control system achieves 
a stable resonance. The optics and the servo loop that controls their positions form 
the core systems of the interferometer; however the many subtleties involved in the 
design and operation of these detectors are outside the scope of this thesis. 

2.2.2 Noise sources in an Interferometer 

In reality, there are many sources of noise which can result in an apparent phase shift 
of the laser light. We define the interferometer strain signal, s, to be the relative 
change in the lengths of the two arms of the interferometer 

s{t) = ^- ^. (2.52) 

Ij 

This signal has two major additive components: (i) a gravitational wave signal h{t) 
and (ii) all other noise sources n{t). The task of gravitational wave data analysts is 
to search for astrophysical signals hidden in this data. The primary goal of the exper- 
imenters engaged in commissioning the LIGO detectors is the reduction of the noise 
appearing in s{t). The noise in interferometers is measured as the amplitude spectral 
density h{f). This is the square root of the power spectral density of the interfer- 
ometer strain in the absence of gravitational wave signals. Figure [5] shows the target 
noise spectral density of the initial LIGO detectors. There are three fundamental 
noise sources that limit the sensitivity of these detectors: 

1. Seismic noise. This is the dominant noise at low frequencies, / ^ 40 Hz. 
Seismic motion of the earth couples through the suspensions of the mirrors and 
causes them to move. To mitigate this, a system of coupled oscillators is used 
to isolate the mirror from the ground motion. 

2. Suspension thermal noise. This noise source limits the sensitivity of the inter- 
ferometer in the range 40 Hz < / < 200 Hz. The steel wire suspending the 
mirror is at room temperature and thermal motion of the particles in the wire 
produce motion of the mirror and change the arm length. 

3. Photon shot noise. At high frequencies, / > 200 Hz, the noise is dominated by 
the shot noise due to the photon counting statistics discussed in the previous 
section. 



For a detailed review of the noise sources present in LIGO's kilometer scale interfer- 
ometers, we refer the reader to [21]. 



2.2.3 Calibration of the Data 

We do not directly record the interferometer strain s{t) but rather the error signal 
v{t) of the feedback loop used to control the differential lengths of the arms. This 
signal, designated LSC-AS_Q in LIGO, contains the gravitational wave signal along 
with other noise. The interferometer strain is reconstructed from the error signal 
in the frequency domain by calibrating v{t) using the response function R{f) of the 
instrument: 

~s{f) = Rimf). (2.53) 

The response function depends on three elements of the feedback control loop shown in 
figure |6] the sensing function C(/); the actuation function A{f ); the digital feedback 
filter D{f)^. 

The sensing function C{f) measures the response of the arm cavities to gravita- 
tional waves. It depends on the light power in the arms, which changes over time as 
the alignment of the mirrors change. The actuation function A{f) encodes the dis- 
tance the mirrors move for the applied voltage at the electromagnets. The dominant 
contribution to this is the pendulum response of the suspended mirrors. The digital 
filter D{f) converts the error signal v{t) into a control signal that is sent as actuation 
to the mirrors to keep the cavities resonant. 

If g{f) is the Fourier transform of the control signal applied to the mirrors, then 
the residual motion of the mirrors is given by 

r{f) = ~s{f) - AUYqU) (2.54) 
as seen in figure [6j The corresponding error signal is 

vU) = C(/)r(/) (2.55) 
and following around the servo control loop we obtain 

~g{f) = DU)vU) = D{f)C{mf). (2.56) 



Substituting equation (2.56) into equation (2.54) and solving for r(/), we obtain 
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where G{f) is the open loop gain of the interferometer, defined by G{f) — A{f)D{f)C{f). 
The error signal is then 

^(/) = C(/)f(/) = 5(/)^^^ (2.58) 

and hence 

«(/) = (2.59) 

The value of the digital filter D{f ) is a known at all times. The actuation function 
can be measured by configuring the interferometer as a simple Michelson, driving a 
mirror and counting the number of fringes that appear at the photodiode for a given 
applied signal. This provides a measure of the displacement of the mirror for a 
given control signal. Since A{f) is due to the pendulum response of the mirror and 
known filters used in the electronics that drive the motion of the mirror, it does not 
change and its value can be established before data taking. A sinusoidal signal of 
known amplitude that sweeps up in frequency is added to the control signal after 
the interferometer is brought into resonance. By comparing the amplitude of this 
calibration sweep in the output of the detector to the known input, the value of the 
open loop gain and (and hence the sensing function) can be determined as a function 
of frequency. The values of the sensing function and open loop gain at the time of 
calibration are denoted Co(/) and Go(/)- 

Although the LIGO detectors have an alignment control system that tries to keep 
the power in the arms constant, the power in the cavity can still change significantly 
over the course of data taking. These fluctuations in power mean that the sensing 
function can change on time scales of order minutes or hours. To measure C(/) 
during data taking sinusoidal signals of known amplitude and frequency /cai added 
to the control signals that drive the mirrors. These calibration signals show up as 
peaks in the spectrum and are called calibration lines. By measuring the amplitude 
of a calibration line over the course of the run compared to the time at which the 
calibration sweep was taken, we may measure the change in the sensing function 

C(f;t)^a(t)Co(f), (2.60) 



where a{t) is the ratio of the calibration line amplitude at time t to the reference 
time. We also allow the digital gain of the feedback loop to vary by a known factor 
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f3{t) so 

D{f;t)=P{t)Do{f). (2.61) 
The response function at any given time, t, becomes 

To analyze the interferometer data we therefore need the error signal, which 
contains the gravitational wave signal, the functions Co(/) and Go{f), which contain 
the reference calibration, and the values of a{t) and P{t), which allow us to properly 
calibrate the data. 



2.3 Gravitational Waves from Binary Inspiral 

Consider a circular binary system comprised of two black holes mi,m2 ~ Mq, sepa- 
rated by a distance a. If a ^ 2GM/c^, where M = mi + m2, then Newtonian gravity 
will provide a reasonably accurate description of the binary dynamics. If we neglect 
higher order multipoles, the gravitational wave field is determined by the quadrupole 
formula [T3] 

2G <PtJ-{t-r) 

where Xjk is the quadrupole moment of the binary, defined by 

X,.^/pW.,..d'. (2.64) 

and is the transverse traceless part of of Ij^. Since the binary can be described 
by Newtonian theory, Kepler's laws are satisfied and the orbital angular velocity is 



GM 

^ = \—r- (2-65) 



a3 



In a Cartesian coordinate system (x, y, z) with origin at the center of mass of the 
binary as shown in figure [7| the mass distribution of the binary is given, in the point 
mass approximation, by 

p{x) = nil — cos VLt)5{ii — Ti sin VLt)5{z)\ 

+ m2 [5{x + r2 cos flt)6{y + r2 sin Qt)6{z)] , 
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where 



ri = a, (2.67) 

mi + 1712 
mi 

ra = — a. 2.68 

mi + m2 

Introduce a second, spherical polar coordinate system labeled (r, t, 0o) related to the 
Cartesian coordinate system by 

et = cos L cos 4>o^x + cos L sin ^o^j) — sin te^, (2.69) 
e^^ = - sin (poBs: + cos 0oey. (2.70) 

To calculate the gravitational radiation hj^ seen by an observer at a position (r, t, 0o) 
relative to the center of mass of the binary we first calculate the quadrupole moment 
of the mass distribution in the frame of the binary. The non-zero components of 



xy 



(2.64) gives 



T 



Xya;. The detailed derivation of X^^, following from equation 

nil [(^{x — ri cos ^t)6{y — ri sin VLt)5{z)] 

+ m2 [5{x + r2 cosfit)5(y + r2 sinnt)5(z)] ci^a; 
(mir^ + m2r2) cos^ VLt 



mi 



m2 



mi + m2 



+ m2 



mi 



mi + m2 



a cos VLt 



Here 



17111712 + m2m\ 
(mi + m2)2 
/ mim2(mi + m2 
V (mi + m2)2 

/ia^ cos^ f2t 

-/xa^ (1 + cos2fit) . 



(2.71) 



a? cos^ fit 
cos^ Vtt 



17111712 



is the reduced mass of the binary and we have used 



5{x - xo)f{x) d^x = /(xo), 



(2.72) 



(2.73) 



6{x) d^x 
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and 



The other components are derived in a similar way to give 

Xyy = ^//a^ (1 - cos 2Q.t) {2.1 A) 

Xxy = Tyx = -fia^ sin 2Qt. (2-75) 
The second time derivative of the quadrupole moment is then 

X,^ = -2/xa^Q^ cos 2Qt (2.76) 

Xyy = 2ixa^^^ cos 2fit (2.77) 



-2fjLa^n^ sm2nt (2.78) 



in the frame of the binary. We transform these to the frame of the observer using the 
standard relations 

_ dxk dxi 

e:=|fe. (2.80) 

to obtain 

Xi, = ixx cos^ i cos^ 00 + iyy cos^ i sw? (j)Q + 2Ixy cos^ i sin 00 cos 00 

= — 2//a^Jl^ cos(2Qi) cos^ i cos 20o — 2//a^Q^ sin(2Qi) cos^ u sin 20o (2-81) 
— — 2//a^Q^ cos^ i cos {2Q:t — 20o) ■ 

Similar transformations give the other components of Xj^ 

X^Q(j)Q = 2/ia^r2^ cos(2r2t) cos20o + -jjio^^^ sin(2r2t) sin20o 

= 2na^VL^ cos (20i - 20o) , (2.82) 
ii^^ = i^^i = -2fjLa^n^ COS t sin (2Qi - 20o) . (2.83) 

Since these are the transverse components of Xij, we can simply remove their trace 
to obtain 2^^: 

7-TT _ _fTT _ f - - ff \ - - ff -T 

= -//a^n^ (l + cos^ l) cos {2nt - 20o) , '^^"^^^ 
= = 2//a2^^2 cos t sin(2Qi - 20o) 
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It is now a simple matter to compute the form of the gravitational radiation using 



the quadrupole formula in equation (2.63) 



hr = = -?^^(l + cos^ cos {2nt - 20o) (2.85) 



= = -^^^ cos . sin {2nt - 20o) • (2.86) 

We can further simplify these equations by using Kepler's third law 

'GMY 

-w) P-**^' 

and defining the gravitational wave frequency f which is twice the orbital frequency 

/ ^ 2 (il) . (2.88, 

If we use the basis vectors and e^^ as the polarization axes of the gravitational 
wave we obtain 

2G 2 

hJt) = hl^ = — — /i(7rGM/)s(l + cos2i)cos(27r/t-20o) (2.89) 
hx{t) = /iJJ = -^K^GMf)! cosisin(27r/t - 20o) (2.90) 

If the vector e^, points from the binary to our gravitational wave detector, the angles l 
and 00 are known as the inclination angle and the orbital phase and r is the luminosity 
distance from the detector to the binary. 

We assume that the binary evolves through a sequence of quasi-stationary circular 
orbits. The orbital energy for a binary with given separation, a, is given by the 
standard Newtonian formula 

IGfiM 

E = —. 2.91 

2 a ^ ^ 



The loss energy due to quadrupolar gravitational radiation is [15 

dE _ G / d'tJ,^d'll^ \_ 32G'M^fx' 
dt ~ 5c5 \ dt'^ dt^ I ~ 5c5 

and so the inspiral rate for circular orbits is given by 

da dE da QAG^ fiM^ 



(2.92) 



dt dt dE 5c^ a^ 



(2.93) 
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The evolution of a as a function of time can therefore be obtained by integrating 

da = ^/iM^ dt (2.94) 

fiM\tc-t) (2.95) 



-^-^^M^] {t,-tY (2.96) 



4 5c5 

and so 

which tells us that the orbit shrinks as orbital energy is lost in the form of gravi- 
tational waves. As the orbit shrinks, the orbital frequency increases and hence the 
gravitational wave frequency and amplitude increase. We call this type of evolution 
a chirp waveform. The evolution of the gravitational wave frequency f{t) can be 
obtained by substituting Kepler's third law, equation (2.87), into equation ( |2.95 ) to 
obtain 



where we have defined rj = fi/M. From this, we may obtain 



(2.97) 



n-i = ^(GM)i^(t.-t). (2.98) 



If we define 6(t) as the dimensionless time variable 



then we obtain 



n-S = (*^]"9(() (2.100) 



which written in terms of the gravitational wave frequency f = Q/n is 



We define the cosine chirp he and the sine chirp hg as 



hcit) = -2[-^] kG'M/(t)]5 cos [2<p{t) - 20o] , (2.102) 



Kit) = ^ ] [rrGMfitp sin [20(t) - 20o] , (2.103) 
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where the orbital phase (j){t) is 

(f){t) = 271 j f{t)dt (2.104) 



and fit) is given by equation (2.101). The + and x waveforms are 

hJt) = -^^(1 + cos^ L)hAt), (2.105) 
hAt) = ^cosihsit). (2.106) 



If the arms of the interferometer form a second Cartesian axis, {x',y',z'), then we 
may define the position of the binary relative to the detector by the spherical polar 
coordinates {r,6,(f). It can be shown that the gravitational waves from the binary 
will produce a strain [23] 

h{t) = F+h+{t) + F^h^it) (2.107) 

at the detector, where the antennae pattern functions and Fx of the detector are 
given by 

F+ = -^(1 + cos^6') cos 2^9 cos 2^ - cos 61 sin 2^9 sin 27/), (2.108) 
Fx = +^(l + cos^6')cos2v9sin2V^-cos6'sin2v9cos2?/'. (2.109) 

The angle if) is the third Euler angle that translates from the detectors frame to the 
radiation frame. The radiation frame is related to the frame of the binary by the 
angles t and 0O) as shown in figure |8] Figure |9] shows the magnitude of the strain 
produced in an interferometer by binary with t = ■?/' = at various positions on the 
sky. It can be seen that the response of the detector is essentially omnidirectional, with 
the maximum sensitivity occurring when the source lies on the 2;-axis of the detector. 
Notice that there are four dead spots in the beam patters where the response of 
the interferometer is zero. These correspond to the locations where the binary is in 
the plane of the interferometer positioned half way between the x and y-axes. We 
will often refer to an optimally oriented binary system. This is a binary located at 
sky position 9 = or tt/2, (i.e. above or below the zenith of the detector) with an 
inclination angle of t = 0. It is so called as this is the position in which the response 
of the detector to the binary is a maximum. 
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2.3.1 Higher Order Corrections to the Quadrupole Waveform 

In the previous section we only considered the lowest order multipole radiation from 
a binary. The goal is to write down a waveform that is sufficiently accurate to 
use matched filtering to search for signals in detector noise. This requires accurate 
knowledge of the phase throughout the LIGO frequency band. In addition to higher 
order multipoles that contribute to the energy loss, there are relativistic corrections 
to the quadrupole formula and effects such as frame dragging and scattering of the 
gravitational wave by the gravitational field of the binary that change the phase 
evolution. Matched filtering is less sensitive to the amplitude evolution, however, so 
we may use the restricted post-Newtonian waveform as the matched filter template. 
The restricted post-Newtonian waveform models the amplitude evolution using the 
quadrupole formula, but includes higher-order v/c corrections to the phase evolution. 
The formula for the orbital phase used in searches for binaries of component mass 
< 3 Mq is given by equation (7) of 



V 



5 



/3715 55 \ ^3 Stt^i 

68 + \ r? 08 64 

V8064 96 7 4 



/ 9 275495 284 875 1855 

+ \ V H V 0* 

V 14 450 688 258 048' 2048'/ 

where 0o and tc are the orbital phase and time at which the binary coalescences and 



9 is defined in equation (2.96). 



2.3.2 The Stationary Phase Approximation 

We will see in chapter [4] that we require the Fourier transforms, hdf) and hs{f), of 
the inspiral waveforms rather than the time domain waveforms given above. In the 
search code, we could compute hc{f) using the Fourier transform of hc{t). This is 
computationally expensive, however, as it requires an additional Fourier Transform 
for each mass pair to be filtered. An alternative method is to use the stationary 
phase approximation [25] to express the chirp waveforms directly in the frequency 
domain [26t [27]. Given a function 

B{t) = A{t) cos 2^{t) (2.111) 

where 

dlnA d(t> , , 
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and 



then the stationary phase approximation to the Fourier transform of B{t) is given by 

1 



B{f) 



m 

A{t) 



exp 



2 27r/t'-20(/)- 



TT 



4 



where t' is the time at which 



Now it is simple to calculate 



exp [-t^if)] 
d(j){t) 
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5 

where we have defined the chirp mass by 

M = fi^Mk 

To obtain the phase function \E'(/) we note that 

ld(j) 
71 dt 



f 



(2.114) 



(2.115) 

(2.116) 

(2.117) 

(2.118) 

(2.119) 
(2.120) 
(2.121) 

(2.122) 
(2.123) 



and we can invert the series in equation (2.110) to write 6 as a function of /. Substi- 
tuting this result and the result for df /dt into the equation for the stationary phase 



approximation, equation (|2.114|), we obtain the form of the inspiral chirps that we 

/GMo' 



will use in matched filtering 



hc{f) 



2GMrr 







5/i 



96M„ 
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(2.124) 
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with the phase evolution given by 

^(/;M,r/)=27r/te-2</)o-7r/4+ ^ 



-5 I / 3715 55 1 3 -|^_^-2 



128// 

15 293 365 27145 3085 ^ , , 

+ I ^ r? H r? x'^ 

^ 508 032 504 ' 72 ' ' 



(2.125) 



where x = (wM fG /c^Y^^ . Notice that in the definition of hc{f) we have absorbed 
the amphtude term 2GMq/ c^r from hj^. This allows us to place at a cannonical 
distance of r = 1 Mpc, as discussed later. Physically the chirp waveform should 
be terminate when the orbital inspiral turns into a headlong plunge, however the 
frequency at which this happens is not known for a pair of comparably massive 
objects. We therefore terminate the waveform at the gravitational wave frequency of 
a test particle in the innermost stable circular orbit of Schwarzschild (ISCQ)[T7] 



which is a reasonable approximation of the termination frequency |2 8]. Since the sine 
chirp is simply the orthogonal waveform to the cosine chirp, we have 

hsU)=^KU)■ (2.127) 



Together equations (2.124), (2.125) and (2.127) give the form of the chirps that we 
will use for the apphcation of matched filtering discussed in chapter |4j 
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Figure 1 : The axes shown in (a) and (b) represent Local Lorentz frames for particles A and C 
respectively. The effect of a gravitational wave in these frame can be described in terms of its effect 
on the vectors ^ and C separating the particles at the origin from particles B and D. 
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Figure 2 : The effect of the two polarizations and /ix of a sinusoidal gravitational wave propa- 
gating through the page on a ring of test particles. As the phase of the gravitational wave changes 
through a complete cycle, the rings are distorted. 
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Figure 3 : The location of the three LIGO interferometers. There are two interferometers at the 
LIGO Hanford Observatory (LHO) in Washington and one interferometer at the LIGO Livingston 
Observatory in Louisiana. 
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Figure 4 : The possible optical configurations of first generation laser interferometers. The inner 
X and y test masses (mirrors) are denoted ITMX and ITMY respectively, and the end x and y test 
masses (mirrors) are denoted ETMX and ETMY respectively. The recycling mirror is denoted by 
RM. Initial LIGO is a power-recycled-Fabry-Perot interferometer, type (c) in this figure. 
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Figure 5 : The fundamental noise sources of LIGO. 
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Figure 6 : The differential more servo control loop. The figure shows the positions of three filters: 
the sensing function C(/), the digital filter D{f) and the actuation function A{f). The input signal 
is s and the measured signal is the error signal v. r is the residual length of the cavity and g is the 
control signal. 




Figure 7 : The parameters of a binary system with rotational axis aligned along the 2;-axis. 
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Chapter 3 

Binary Black Hole MACHOs 

One of the most interesting current problems in astrophysics is that of dark matter. 
Dark matter is so called because it has eluded detection through its emission or 
absorption of electromagnetic radiation. Our knowledge of its existence comes from 
its gravitational interaction with luminous matter in the universe. There have been 
several ideas proposed to explain the nature of dark matter; chief among these are 
weakly interacting massive particles (WIMPs) and massive astrophysical compact halo 
objects (MACHOs) |29] . WIMPs, supersymmetric particles produced as a relic of the 
big bang, are outside the scope of this thesi^ No compelling reason exists to think 
that WIMPs will produce significant gravitational waves. In this chapter, we review 
the evidence for dark matter in the form of MACHOs in the Galactic halo. The 
nature of MACHOs is unknown; we review a proposal that suggests that if MACHOs 
are primordial black holes (PBHs) formed in the early universe, then some of the 
PBH MACHOs may be in binary systems P^. Searching for gravitational waves from 
the inspiral and coalescence of these binary black hole MACHOs (BBHMACHOs) is 
the motivation for this thesis. 

3.1 Dark Matter In The Galactic Halo 

Dark matter is detected by its gravitational interaction with luminous matter. Strong 
evidence for the presence of dark matter in the universe comes from the study of 
galactic rotation curves: measurements of the velocities of luminous matter in the 
^We refer to |30j for a review of the nature of dark matter. 
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disks of spiral galaxies as a function of galactic radius. Consider a simple rotational 
model for the disk of a spiral galaxy. Consider a star with mass orbiting at radius 
r outside the disk of the galaxy. Newtonian dynamics tells us that if the mass inside 
radius r is rUg then 

GmglTLs _ ITLsVl 

2 ~ \'^-^) 

where Vg is the velocity of the star and G is the gravitational constant. Let us 
suppose that as we increase r, the change in the irig is negligible, which is a reasonable 
assumption towards the edge of the disk of a typical spiral galaxy. We can see from 



equation (3.1) that we would expect the velocity of stars at the edge of the galactic 
disk to fall off as 

1 , , 

Vs oc (3.2) 



when r ^ R, where R is the radius containing most of the disk matter. Galactic 
rotation curves, determined using the Doppler shift of the 21 cm hydrogen line, have 
been measured for several galaxies[3T]. It is found that the rotation curves do not fall 
off as expected. Instead the rotational velocities of galactic matter are measured to 



be constant out to the edge of the visible matter in the disk, as shown in figure 10 



This surprising result suggests that 80%-90% of the matter in spiral galaxies is in the 
form of dark matter stretching out at least as far as the visible light. 

A typical argument to understand the formation of galactic disks from baryonic 
matter considers an initially spherical distribution of baryonic matter rotating with 
some angular momentum, L. Over time the matter will lose energy through inelastic 
collisions. Since the angular momentum of the system is conserved, the initial distri- 
bution must collapse to a rotating disk. On the other hand, if the initially spherical 
distribution is composed of dark matter instead of baryons, the collisions will be elas- 
tic because the dark matter is weakly interacting. As a result of this, dark matter 
initially distributed in an isotropic sphere will maintain this distribution over time. 
Since we do not expect a spherical dark matter halo to collapse to a disk, the sim- 
plest possible assumption is that the dark halo is a spherical, isothermal distribution 
of dark matter. This suggests that dark matter will be distributed in an extended 
halo encompassing the luminous matter of a galaxy. If we assume that the density of 
the dark matter is p(r) then the mass within a thin shell of a spherical halo is 

dM{r) = 47rrV(r) dr, (3.3) 



where dr is the thickness of the shell. Using Newtonian dynamics, the velocity w of a 
particle of mass m at radius r is 

GM{r)m mv"^ 

^ (3.4) 

r 

The galactic rotation curves tell us that the velocity is independent of the radius, so 

2 

V r 



M{r) = — . (3.5) 



Differentiating this with respect to r and substituting the result into equation (3.3), 
we obtain 

dM{r) ^ ^ . . ,^ ^, 

= - = 47rrV(r) (3.6) 

which gives 

"<■'■) = <^-^' 
If we assume that the dark and visible matter are in thermal equilibrium, we may 
use the measured rotational velocity of local stars about the galactic center as the 
velocity of the dark matter. 

We can easily estimate the density of dark matter in the neighborhood of the 
Earth p(r£;) as follows. The earth is approximately 8 kpc from the galactic center 
and the rotational velocity of objects at this radius is f ~ 200 kms~^. Using these 



values in equation (3.7), we find 

p(rs) = 7.6 X 10-2^gcm-3. (3.8) 

More sophisticated modeling of the Galaxy [32], suggests that the local halo density 
is 

p(rs) = 9.2+3-8 ^ i0"25gcjn-3 (3.9) 
or approximately 0.01 Mq pc"^. 



Equation (3.7) applies only at intermediate radial distances. The data at small 
r is consistent with the dark matter having a constant core density pc within a core 
radius a[H3]. The halo density then becomes 
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The values of pc and a are obtained by fitting measured galactic rotation curves to 



equation (3.10) using data near the galactic center. There is, in fact, no evidence to 
suggest that halos are exactly spherical. In fact the halo density may be flattened |33]. 
For a flattened halo a model of the dark matter density becomes 

where R and z are galactocentric cylindrical coordinates and g is a parameter that 
describes the flattening of the halo. At present there is no measurement of the extent 
of galactic halos beyond the luminous matter. For the Milky Way it is thought that 
the halo extends out to a radius of ~ 50kpc, although it is possible that it extends 
all the way out to the Andromeda galaxy at ~ 700 kpc. 



3.2 MACHOs in the Galactic Halo 

Galactic rotation curves provide strong evidence that spiral galaxies such as the Milky 
Way are surrounded by a large quantity of dark matter, but tell us nothing about the 
nature of this dark matter. A variety of candidates have been proposed to explain the 
nature of dark matter. These generally fall into two classes. The flrst class consists 
of elementary particles such as axions|Mj or weakly interacting massive particles 
(WIMPs) • Such dark matter candidates are outside the scope of this thesis. Active 
searches for WIMPs and axions are underway and we refer to [30] for a review of the 
particle physics dark matter candidates. The second class of dark matter candidates 
are known as massive astrophysical compact halo objects or MACHOs. MACHOs 
are objects such as brown dwarfs (stars with insufficient mass to burn hydrogen), red 
dwarfs (stars with just enough mass to induce nuclear fusion), white dwarfs (remnants 
of 1-8 Mq stars) or black holes located in the halos of galaxies. Optical and infrared 
observations in the early 1990's were not sensitive enough to constrain the fraction of 
the halo in MACHOs [3^] and the method of gravitational lensing was suggested as a 
method for detecting halo dark matter in the form of MACHOs|37]. 
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3.2.1 Gravitational Lensing of Light 

Gravitational lensing is caused by the bending of light around a massive object. As- 
sume that a MACHO produces a spherically symmetric gravitational field; the geom- 
etry of spacetime around the MACHO satisfies the Schwarzschild solution. Consider 



the scattering of light by a MACHO shown in figure 11 , where b the impact parameter 
of the light, the minimum distance of the photon to the MACHO. Recall that the 
lightlike orbits of Schwarzschild spacetime satisfy pTTj 

d'u 3GM 2 

where u = b/r and M is the mass of the MACHO. If R is the size of the MACHO, 
then 

3GMu^ 3GMR 



if 6 ^ -R. we can solve equation (3.12) perturbatively in the small parameter e = R/b 
as follows. Write 

u = Uq + eui + • • • (3-14) 



and substitute into equation (3.12) to get 

Un + eu'l-l h Mo + ewi H = (mq + ewi H )^ (3.15) 

c^R 

where prime denotes differentiation with respect to 0. At leading order, 

u'^ + uo = (3.16) 

has solution 

Uo = Asm{(f) + (j)o) (3.17) 
where A and (pQ are constants. We are free to choose any value for (pQ as it simply 



chooses an orientation for the axes in figure 11, so let 0o = 0. Since (f) = tt/2 gives 
the distance of closed approach, we find A = b/r^i^i = 1- Now ui satisfies 

3GM . 

u^+ui= sm 

3GM , ^^-^^^ 
^1 — cos 20) . 



2c^R 

Inspection suggests a solution of the form 



3GM 

•"1 = o 9 n + cos 20. (3.19) 
Ic^R 
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Substituting this into equation (3.18) we find that 

, 3GM , 3GM 3GM 

- 4a cos 2^ + — + a cos 2^ = — - — cos 2^ (3.20) 

and so a = GM /2c^R. The solution for m, up to first order, is therefore 

GMe , , , , 

u = sin0+— -(3 + COS20). (3.21) 

As r ^ oo, u ^ and —5/2, so 

= sin (^-- j + — (3 + cos(-5)) (3.22) 

For small 5, sin(— 5/2) —S/2 and cos(— 5) ~ 1, so 

5 2GM 

0«-- + ^ (3.23) 

. 4GM 

6 ^ — 3.24 



The total defiection of the light is therefore 

AGM 



6 = — . (3.25) 



Suppose a MACHO lens is at a distance Dsl from a source star and an observer is at 



a distance Z^l from the MACHO as shown in figure 12 Then the ray of light from 
the source that encounters the MACHO with critical impact parameter rg will reach 
the observer. Simple geometry, using the small angle approximations, shows that 

rp rp AGM 
5 = ^s + ^o = 7f + 7^ = ^ (3.26) 

therefore the observer sees the lens light when the ray is at the Einstein radius, te, 
given by 

/4GM JslDl ' ,„ 



If the source, MACHO and observer are coUinear, as shown in figure 12 the observer 
sees a bright ring of radius te around the MACHO. The angular radius of this ring 
is the Einstein angle, 

Iagm Id^i, 
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In the realistic case of slight misalignment, then the lensed star will appear as two 
small arcs. Consider a MACHO of mass 0.5 Mq at a distance of = 25 kpc lensing 
a star in the Large Magellanic Cloud (LMC) at a distance of 50 kpc. Then 



2GM 



9^=^ =^ ~ 10-^° ^ 2" X 10-^ (3.29) 

too small to be resolved by optical telescopes. Fortunately the lensing produces an 
apparent amplification of the source star by a factor [38] 

A = ""l^ , (3.30) 

where v = 13/6-^, and (3 is the angle between the observer-lens and observer-star lines. 
Since objects in the halo are in motion, 

m = yM)^+^, (3.31) 

where v± is the transverse velocity of the lens relative to the line of sight, f3^i^ is 
the closest approach of the lens to the source, and t is the time to the point of 
closest approach [37t 129] . Searches for the amplification of stars caused by gravitational 
lensing of 6'e ~ micro arc seconds are referred to as gravitational microlensing surveys. 
Such surveys measure magnification of the star and the duration of the microlensing 
event. Unfortunately it is not possible to determine the size of the the lens from these 
measurements. 

3.2.2 Gravitational Microlensing Surveys 

Several research groups are engaged in the search for microlensing events from dark 
matter |39l HQ]. By monitoring a large population of well resolved background stars 
such as the LMC, constraints can be placed on the MACHO content of the halo. The 
MACHO project has conducted a 5.7 year survey monitoring 11.9 million stars in the 
LMC to search for microlensing events (39j using an automated search algorithm to 
monitor the light curves of LMC stars. Optimal filtering is used to search for light 



curves with the characteristic shape given by equation (3.30). 

Since the effect of microlensing is achromatic, light curves are monitored in two 
different frequency bands to reduce the potential background sources which may 
falsely contribute to the microlensing rate. Background events include variable stars 
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in the LMC (known as bumpers |41]). which can usually be rejected as the fit of the 
light curves to true microlensing curves is poor. Supernovae occurring behind the 
LMC are the most difficult to cut from the analysis. The MACHO project reported 
28 candidate microlensing events in the 5.7 year survey of which 10 were thought to 
be supernovae behind the LMC and 2-4 were expected from lensing by known stellar 
populations. They report an excess of 13-17 microlensing events, depending on the 
selection criteria used. 

The optical depth, r, is the probability that a given source star is amplified by a 
factor A > 1.34 [37]. This is just the probability that the source lies on the sky within 
a disk of radius Be around a microlensing object and is given bv|42j 



^ = — / P{l)^-r^dl, (3.32) 
c ./o 



where L = Dsl + Dl is the observer-star distance and / = is the observer-lens 



distance. For the spherical halo given in equation (3.10) with density 



p2 _|_ 2 

p(r) = 0.0079^^ ^ Mq pc-^ (3.33) 

-|- a'' 

where i?o = 8.5 kpc is the Galactocentric radius of the Sun and a galactic core radius 
of a = 5 kpc, the predicted optical depth towards the LMC (assumed to be at 50 kpc) 
is [12] 

tlmc = 4.7 X 10-^. (3.34) 

The optical depth towards the LMC measured by the MACHO project microlensing 
surveys is 

TLMC = l.2t^oi X 10-^ (3.35) 

This suggests that the fraction of the halo in MACHOs is less that 100%, but does 
not exclude a MACHO halo. 

The number of observed MACHO events and the time scales of the light curves 
can be compared with various halo models. The MACHO project has performed 
a maximum-likelihood analysis in which the halo MACHO fraction / and MACHO 
mass m are free parameters. For the standard spherical halo, they find the most likely 
values are / = 20% and m = 0.45 Mq. The 95% confidence interval of on the MACHO 
halo faction is / = 8%-50% and the 95% confidence interval of on the MACHO mass 
is 0.15-0.9 M0. The total mass in MACHOs out to 50 kpc is found to be 9^3 x 
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10^° Mq, independent of the halo model[39j. The EROS collaboration has recently 
published results of a search for microlensing events towards the Small Magellanic 
Cloud (SMC)[40J. The EROS result further constrains the MACHO fraction of a 
standard halo in the mass range of interest to less than 25%; they do not exclude a 
MACHO component of the halo, however. 

3.3 Gravitational Waves from Binary Black Hole MACHOs 

Since the microlensing surveys have shown that ~ 20% of the halo dark matter may 
be in the form of ~ 0.5 Mq MACHOs, it is natural to ask what the MACHOs may 
be. As we mentioned above, it has been proposed that MACHOs could be baryonic 
matter in the form of brown dwarfs, objects lighter than ~ 0.1 Mq that do not have 
sufficient mass to sustain fusion, however, this is inconsistent with the observed masses 
of MACHOs. The fraction of the halo in red dwarfs, the faintest hydrogen burning 
stars with masses greater than ~ 0.1 Mq, can be constrained using the Hubble Space 
Telescope. Hubble observations may also be used to constrain the fraction of the halo 
in brown dwarfs. The results indicate that brown dwarfs make up less than ~ 3% and 
red dwarfs less than ~ 1% of the halo[l3l HI]. A third possible candidate for baryonic 
MACHOs is a population of ancient white dwarfs in the halo. White dwarfs are the 
remnants of stars of mass 1-8 Mq and have masses of ~ 0.6 Mq. Although they seem 
to be natural candidates for MACHOs, proper motion searches for halo white dwarfs 
have been conducted and no candidates have been found|l5l SHI H?]. Creeze et al. 
combined the results of previous surveys to find that 4% (95% confidence level) of 
the halo is in the form of white dwarfs |17]. 

It is possible that there is an over dense clump of MACHOs in the direction of 
the LMCjlH], the lenses are located in the LMC itself |19] or the lenses are in the 
disk of the galaxy [50]. If the MACHOs detected by microlensing are truly in the 
halo, however, it is possible that MACHOs are non-baryonic matter such as black 
holes [511 [13] • Black holes of mass ~ 0.5 could not have formed as a product 
of stellar evolution and so they must have been formed in the early universe [521 153] . 
Several mechanisms have been proposed to form primordial black holes with the 
masses consistent with the MACHO observations. These include multiple scalar fields 
during inflation[51], chaotic inflation[55j or reduction of the speed of sound during 
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the QCD phase transition [56]. We do not consider these formation mechanisms in 
detail here; it is sufficient for our purposes that PBHs with masses consistent with 
microlensing observations can form. If the MACHOs are primordial black holes then 
there must be a large number of them in the halo. The total mass in MACHOs out 
to 50 kpc is 9 X 10^° M0, as measured by microlensing surveys. If these are 0.5 
PBHs then there will be at least ~ 1.8 x 10^^ PBHs in the halo. With such a large 
number of PBHs in the halo it is natural to assume that some of these may be in 
binary systems. 

Nakamura et al. [13] considered PBHs formed when the scale factor of the universe 
R, normalized to unity at the time of matter-radiation equality, is 

Rf = ^/GM^c^L,^ = 1.2 X 10-« ' ^3.36) 

where Lgq is the Hubble horizon scale at the time of matter-radiation equality, Vt is 
the fraction of the closure density in PBHs and h is the Hubble parameter in units of 
100 km s^^. The age and temperature of the universe at this epoch are ~ 10~^ seconds 
and ~ 1 GeV, respectively. By considering a pair of black holes that have decoupled 
from the expansion of the universe to form a bound system interacting with a third 
black hole, which gives the pair angular momentum to form a binary, they showed 
that the distribution of the semi-major axis, a, and eccentricity, e of a population of 
binary black hole MACHOs is 

f{a,e)dade= — ^dade (3.37) 

2x5(1 — 6^)2 

where x is the mean separation of the black hole MACHOs at the time of matter- 
radiation equality, given by 

X = 1.1 X 10^^ f — J [nh'^) 5 cm, (3.38) 

The coalescence time of a binary by the emission of gravitational waves is approxi- 
mately given by [57] 



t = to(^) (l-e')^ (3.39) 
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where to = 10^° years and 

ao = 2 X 10^^ [H-] cm (3.40) 



M \3 



Mr. 







is the semimajor axis of a binary with circular orbit which coalesces in time to- 



Integrating equation (3.37) for fixed t using equation (3.39), Nakamura et a/. [13 



found the probability distribution ft{t) for the coalescence time is 



_3_ 3' 



^max / V ^max 



dt 

J. (3.41) 



where tmax = tQ{x/ao)^. The number of coalescing binaries with t ~ to is then 
~ 5 X 10® for Qh"^ = 0.1, so the event rate of coalescing binaries is ~ 5 x 10^^ events 
per year per galaxy. loka et al. |58] performed more detailed studies of binary black 
hole MACHO formation in the early universe and found that, within a 50% error, the 
distribution function and the rate of coalescence given in [I3| agree with numerical 
simulations. The event rate of coalescing binary black hole MACHOs is therefore 

i?BBHMACHO = 5 X lO'^ X 2^^ ^1'^ galaxy-\ (3.42) 

This rate is significantly higher than the coalescence rate of binary neutron stars, 
which isjSU] 

i?BNS = 8.3 X 10-^ X 2^^ yr"^ galaxy'^ (3.43) 

It must be emphasized that several neutron star binaries have been observed, but 
there are no observations of black hole MACHO binaries. 

The distance to which we can detect a binary inspiral is usually expressed in terms 
of the characteristic strain, h^har of the binary. This represents the intrinsic amplitude 
of the signal at some frequency times the square-root of the number of cycles over 
which the signal is observed at that frequency 



KtUf) = \fHf)\ - (3.44) 
For an inspiral signal this is given by[16j 



-1 



W/) = 4xlO-^-j ^— j , (3.45) 

where r is the distance to the binary and Ai is the chirp mass. For comparison with 
signal strength, the detector sensitivity is better expressed in terms of the root mean 
square (RMS) dimensionless strain per logarithmic frequency interval 



where s is the detector strain output in the absence of a gravitational wave signal and 
Sn{f) is the power spectral density of s. If the value of /ichar > (a few) x /irms, then 
the binary will be detectable. Figure [13] shows the characteristic strain of a binary 
consisting of two 0.5 Mq black holes at 10 Mpc compared to the RMS noise for initial 
LIGO. It can be seen that the inspiral signal lies significantly above the noise, so 
these MACHO binaries could be excellent source for initial LIGO. Nakamura et a/.|13j 
showed that the rate of MACHO binaries could be as high as 3 yr~^ at a distance of 
15 Mpc, under their model assumptions. 

3.4 Binary Black Hole MACHO Population Model 

The goal of this thesis is to search for gravitational waves from binary black hole 
MACHOs described in the previous section. In the absence of a detection, however, 
we wish to place an upper limit on the rate of binary black hole MACHO inspirals 
in the galaxy. We can then compare the predicted rate with that determined by ex- 
periment. We will see later that in order to determine an upper limit on the rate, we 
need to measure the efficiency e of our search to binary black hole MACHOs in the 
galactic halo. We do this using a Monte Carlo simulation which generates a popula- 
tion of binary black hole MACHOs according to a given probability density function 
(PDF) of the binary black hole MACHO parameters. Using the set of parameters 
generated by sampling the PDF, we can simulate the corresponding inspiral wave- 
forms on a computer. We then digitally add the simulated waveforms to the output 
of the gravitational wave detector. By analyzing the interferometer data containing 
the simulated signals, we can determine how many events from the known source 
distribution we find. The efficiency of the search is then simply 
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Recall that an inspiral waveform is described by the following nine parameters: 
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To simulate a population of BBHMACHOs in the halo we need to generate a list of 
these parameters that correctly samples their distributions. 

We first address the generation of inspiral end time, t^. The nature of the noise 
in the interferometers changes with time, as does the orientation of the detectors 
with respect to the galaxy as the earth rotates about its axis over a sidereal day. To 
sample the changing nature of the detector output, the Monte Carlo population that 
we generate contains many inspiral signals with end times distributed over the course 
of the science run. We generate values of tc at fixed intervals starting from a specified 
time to- The fixed interval is chosen to be 2048 + 7r 2051.141592653 sec. This allows 
us to inject a significant number of signals over the course of the two month run with 
the signals far enough apart that they do not dominate the detector output. The 
interval is chosen to be non-integer to avoid any possible periodic behavior associated 
with data segmentation. The start time for the Monte Carlo, to? is chosen from a 
uniform random distribution in the range tstart — 2630/7r < < tstart; where tgtart 
is the time at which the science run begins. We stop generating inspiral parameters 
when tc > tend the time at which the second science run ends. For each generated 
end time, tc we generate the other inspiral waveform parameters. 

We obtain the distribution of the mass parameters {1711,7712) from the microlensing 



observations of MACHOs in the galactic halo, described in section 3.2.1 , which suggest 
that the most likely MACHO mass is between 0.15 and 0.9 Mq. In the absence of 
further information on the mass distribution we simply draw each component mass 
from a uniform distribution in this range. We increase the range slightly to better 



48 



measure the performance of our search at the edge of the parameter space. We also 
note that the search for binary neutron stars covers the mass range 1 to 3 M0, so we 
continue the BBHMACHO search up to 1 Mq rather than terminating it at 0.9 Mq. 
We therefore generate each BBHMACHO mass parameter, mi or m2, from a uniform 
distribution of masses between 0.1 and 1.0 Mq. 

The angles l and 0o are generated randomly to reflect a uniform distribution in 
solid angle; cost is uniform between —1 and 1 and 0o is uniform between and 2n. 
The polarization angle ip is also generated from a uniform distribution between and 
27r. 

To generate the spatial distribution of BBHMACHOs, we assume that the dis- 
tribution in galactocentric cylindrical coordinates, {R,6,z), follows the halo density 



given by equation (3.11), that is, 

p(r) oc TTi-^ (3.48) 

where a is the halo core radius and q is the halo flattening parameter. We can see 
that this distribution is independent of the angle ^, so we generate Q from a uniform 
distribution between and 2i\ . If we make the coordinate change zjq —>■ z, we 
may obtain a probability density function (PDF) for the spatial distribution of the 
MACHOs given by 

f(R,z) RdRdz oc — ^ 7;RdRdz. (3.49) 

a'' + R^ + z^ 

We wish to randomly sample this PDF to obtain the spatial distribution of the BBH- 
MACHOs. Once we have obtained a value of the new coordinate z, we simply scale 
by q to obtain the original value of z. Recall that for a probability density function 
/(x) the cumulative distribution FiJC) given by 

F{X) = j f{x) dx (3.50) 

with F{X) G [0, 1] for all f{x). If we generate a value of u from a uniform distribution 
between [0,1] and solve 

f{x')dx' = u (3.51) 



for X then we will uniformly sample the probability distribution given by f{x). Notice, 



however, that PDF in equation (3.49) is a function of the two random variables R and 
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z, rather than a single variable x as in equation (3.51). Let us make the coordinate 
change 



R = r cos ip, 
z = r sin (p 



(3.52) 



and so 



tancp 



z 

R^ + z\ 



(3.53) 



Equation (3.49) becomes 
1 



/C 





r cos 09 



o „o o RdRdz = K, 
+ R^ + z^ J J a'^ + r 

1 pr 

I Jo 



/C [sin 'p'^_^ 



g2 _|_ J.2 



dr 



(3.54) 



r — a arctan 



where r„ 



50 kpc is the extent of the halo and /C is a constant that normalizes the 



PDF to unity. We can see immediately from equation ( |3.54[ ) that sin ip is uniformly 

(3.55) 



distributed between —1 and 1. Now consider the PDF for r given by 

'r 



/(r) = /C 



r — a arctan 



with normalization constant 



-^max ^ circtjcLll 



R„ 



(3.56) 



To sample the distribution for r, generate a random variable u uniform between 
and 1 and find the root of 



r — a arctan 



u 



R„ 



a arctan 



0. 



(3.57) 



We can see that the value of r that solves equation (3.57) must lie between and 



i?max and that the left hand side is a monotonically increasing function of r. We may 
therefore use a simple bisection to solve for the value of r. The values of r and ip are 



easily inverted for R and z using equation (3.52) 
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This method was implemented in lalapps_minj and figure 14 shows a histogram of 
the first mass parameter generated the by the Monte Carlo code. It can be seen that 



this is uniform between 0.1 and 1.0 Mq, as expected. Figure 15 shows the spatial 
distribution of BBHMACHO binaries for a spherical, g = 1, halo that extends to 
-Rmax = 50kpc with a core radius of a = 5 kpc. Since the software that simulates 
inspiral waveforms expects the position of the inspiral to be specified in equatorial 
coordinates, the population Monte Carlo code also generates the coordinates of the 
inspiral as longitude, latitude and distance from the center of the earth, as shown in 
figure 16 We will return to the use of population Monte Carlos in chapter [7j 
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Figure 10 : The results of observations of the hydrogen 21 cm line of the spiral galaxy NGC 3198 
show that the rotation curve is flat out to the last measured point at 30 kpc[60 . This implies a 
large discrepancy between the observed rotation curve and that predicted from light observations. 
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Figure 11 : A photon will be scattered by the curvature of spacetime caused to the gravitational 
field of a MACHO. If the closest approach of the photon is b, it will be deflected by an angle 
6 = AGM/(?b. 




Figure 12 : The geometry of microlensing of light from a source star by a MACHO showing the 
definition of the Einstein radius, Te- 



53 



10-i« - 



Co 



10-20 



10-2> - 



10-22 



10-23 



rrm — — i i mini 1 — i mini 1 — i i i 1 1 ii 




LmL 



J I 



J 



J 



l3 



10 



100 



1000 



10* 



/(Hz) 



Figure 13 : The sensitivity of LIGO to a binary black hole MACHO inspiral can be considered in 
terms of comparison between the characteristic strain /ichar of the inspiral with the RMS noise curve 
of the detector. It can be seen that if binary black hole MACHOs exist, they could be an excellent 
source for LIGO. 
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Figure 14 : The BBHMACHO population Monte Carlo code is used to simulate a distribution of 
209048 coalescing binaries and a histogram is made of the first component mass, mi to confirm that 
it is uniformly distributed over the expected range. Similar tests are performed for the second mass 
parameter, m2, the galactocentric longitude, 6, the inclination angle, i, the polarization angle, ip, 
and the coalescence phase, 4>c- 




Figure 15 : The spatial distribution of 5000 simulated BBHMACHO binaries in a spherical q — 1 
Galactic halo of size -Rmax — 50kpc with a core radius a = 8.5 kpc shown in galactocentric coordi- 
nates. Each point in the figure corresponds to a simulated binary black hole MACHO injection. 
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Figure 16 : The spatial distribution of 5000 simulated BBHMACHO binaries in a spherical, q — 0, 
Galactic halo of size i?max — 50kpc with a core radius a — 8.5 kpc shown in equatorial coordinates. 
Each point in the figure corresponds to a simulated binary black hole MACHO injection. The color 
of the point shows the distance from the center of the earth to the binary. Note the dense clump of 
binaries in the southern hemisphere, towards the center of the Galaxy. 
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Chapter 4 

Binary Inspiral Search Algorithms 



Using equation (2.105 )-(2. 109), we may write the gravitational wave strain induced 
in the interferometer as 

where 



V 



cos (20(t) - 9) 



A{t) 



2Gfx 



nGMfit)] 



(4.1) 
(4.2) 



and V is the effective distance, given by 

V 



^Fl{l + C0S2 6)2/4 + F2 C0S2 



The phase angle 9 is 



tan 6' 



Fx 2 cos L 
F+(l + cos2 i) 



(4.3) 



(4.4) 



and 0(t) is given by equation (2.110). In this chapter we address the problem of 
finding such a signal hidden in detector noise. The detection of signals of known 
form in noise is a classic problem of signal processing pT] and has been studied in the 



context of binary inspiral in [621 163]. This material is reviewed in section 4.1 The 
particular implementation used to extract inspiral signals from interferometer data 



in a computationally efficient manner is presented in section 4.3 



4.1 Detection of Gravitational Waves in Interferometer Noise 



Our goal is to determine if the (calibrated) output of the interferometer s{t) contains 
a gravitational wave in the presence of the detector noise described in section 2.2.2 
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When the interferometer is operating properly 

{n(t) + hit) signal present, 
(4.5) 
n{t) signal absent. 

The instrumental noise n{t) arises from naturally occurring random processes de- 
scribed mathematically by a probability distribution function. The optimum receiver 
for the signal h{t) takes as input the interferometer data and returns as its output 
the conditional probability P{h\s) that the signal h{t) is present given the data s{t). 
The conditional probability that the signal is not present, given the data is then 
P(0|s) = 1 — P{h\s). The probabilities P{h\s) and -P(0|s) are a posteriori prob- 
abilities. They are the result of an experiment to search for the signal h{t). The 
probability that the signal is present before we conduct the experiment is the a priori 
probability P{h). Similarly, P(0) = 1 — P{h) is the a priori probability that the 
signal is absent. 

The construction of the optimal receiver depends on the following elementary 
probability theory. The probability that two events A and B occur is given by 

P(A B) = P{A)P{B\A) = P{B)P{A\B), (4.6) 

allowing us to relate the two conditional probabilities by 

, P{A,B) PiA)P{B\A) , , 

If instead of a single event, A, suppose we have a complete set of mutually exclusive 
events Ai,A2,. ■ ■ ,Ak- By mutually exclusive we mean that two or more of these 
events cannot occur simultaneously and by complete we mean that one of them must 
occur. Now suppose B is an event that can occur only if one of the Ak occurs. Then 
the probability that B occurs is given by 

K 



P{B) = Y,P{A,)P{B\A,). (4. 



k=l 



Equation (4.8) is called the total probability formula. Now let us suppose that B is 



the result of an experiment and we want to know the probability that it was event 



Ak that allowed B to happen. This can be obtained by substituting equation (4.8) 
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into equation (4.7) to get 



P{Ak)P{B\A, 
P{B) 



P{Ak)P{B\Au) 
El,PiA,)P{B\A, 



(4.9) 



Equation (4.9) is Bayes' theorem. The probability P{Ak) is the a priori probabihty 
of event Ak occurring and P{Ak\B) is the a posteriori probabihty of Ak occurring 
given that the outcome of our experiment B occurred. The conditional probability 
P{B\Ak) is called the likelihood. 

Now suppose that set {^4^} contains only to two members: "the signal is present" 
and "the signal is absent". The a priori probabilities of these events are P{h) and 
P(0), as discussed earlier. We consider B to be the output of the interferometer for 
a particular experiment. We can use Bayes' theorem to compute the a posteriori 
probability that the signal is present, given the output of the detector: 

P{h)P{s\h) 



P{h\s) 



Pis] 



(4.10) 



where P{s) is the a priori probability of obtaining the detector output and P{s\h) 
is the likelihood function. P{s\h) is the probability of obtaining the detector output 
given that the signal is present in the data. The probability of obtaining the detector 
output is given by 

P{s) = P{h)P{s\h) + P(0)P(s|0) (4.11) 
since the signal is either present or not present. Substituting equation ( 4.11[ ) into 



(4.10), we write 



P{h\s) 



P{h)P{s\h) 



P{h)P{s\h) + P{0)P{s\0)' 



(4.12) 



Dividing the numerator and denominator on the right hand side of equation (4.12) 
by P{h)P{s\0) we obtain 

P{s\h)/Pis\0) 



P{h\s 

Define the likelihood ratio 



[Pis\h)/Pism + [Pio)/p{h)y 

P{s\h) 



A 



so that equation (4.13) becomes 

P{h\s) 



P{s\0) 



A 



A + [p{o)/p{h)y 



(4.13) 



(4.14) 



(4.15) 
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Similarly, we find that the probability that the signal is absent is given by 

P{0\s) = 1 - P{h\s) = ^^"^/f'f^j,^^ . (4.16) 
^ ' ^ ^ ' ^ A+[P{0)/P{h)] ^ ' 



Using equations (4.15) and (4.16), we find that the ratio of the a posteriori probabil- 
ities is 

^ = A^. (4.17) 
P(0|s) P(0) ^ ^ 

We now construct a decision rule for present or absence of the signal. If P{h\s) 

is large (close to unity) then it is reasonable to conclude that the signal is present. 

Conversely, if P{h\s) is small (close to zero) then we may conclude that the signal 

is absent. Therefore we may set a threshold P^, on this posterior probability as our 

decision rule is 

P{h\s) > P=K decide the signal is present, (4-18) 
P{h\s) < P^ decide the signal is not present. (4.19) 

Given this decision rule there are two erroneous outcomes. If P{h\s) > P* and the 
signal is not present, we call this a false alarm; our decision that the signal is present 
was incorrect. Conversely, if P{h\s) < P^, and the signal is present, we have made a 
false dismissal. Each possible outcome has an associated probability 

F probability that we have a false alarm (4.20) 

F' = 1 — D probability that we have a false dismissal, (4-21) 

where D is the probability of a correct detection. To construct the posterior prob- 
ability, P{h\s) we need the unknown a priori probabilities, P{h) and P(0). We see 



from equation (4.15), however, that P{h\s) is a monotonically increasing function of 
the likelihood. The ratio of the a priori probabilities, P{h)/P{0), is a constant that 
does not involve the result of our experiment. Therefore we can define the output 
of our optimum receiver to be the device which, given the input data s{t), returns 
the likelihood ratio A. For the receiver to be optimal in the Neyman- Pearson sense 
the detection probability should be maximized for a given false alarm rate, F. Rule 



(4.18)-(4.19) is optimal in the Neyman-Pearson sense. 

We now consider the construction of A for the interferometer data s{t) and the 
gravitational wave signal h{t). Assume that the noise is stationary and Gaussian with 
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zero mean value 



{n{t)) = 



(4.22) 



where angle brackets denote averaging over different ensembles of the noise. The (one 
sided) power spectral density S'„(|/|) of the noise is defined by 

{n{f)n{f')) = ls^im{f-f) (4.23) 

where n{f) is the Fourier transform of n{t). We wish to compute the quantity 

A = (4.24) 

P(s|0)' ^ ^ 

however since the probabilities P{s\h) and P{s\0) are usually zero, so in calculating 
the likelihood ratio, we must get rid of the indeterminacy by writing 

P{s\h) p{s\h)ds p{s\h) 



A 



(4.25) 



P(s|0) p{s\0)ds p{s\0) 

Instead of using the zero probabilities where P{s\h) and P(s|0), we use the corre- 
sponding probability densities p{s\h) and p{s\0). The probability density of obtaining 
a particular instantiation of detector noise is[62j 

I, , 



p{n) = K, exp 



\n\n] 



(4.26) 



where /C is a normalization constant and the inner product (-I-) is given by 



(a I h) 



df 



Sum 



(4.27) 



The probability density of obtaining the interferometer output, s{t), in the absence 
of signal, i.e. s{t) = n{t), is therefore 



p{s\0) = p{s) = /Cexp 



[S s 



(4.28) 



The probability density of obtaining s{t) in the presence of a signal, i.e. when s{t) 
n{t) + h{t), is given by 



p{s\h) = p{s — h) = /Cexp 



--{s-h\s-h) 



(4.29) 



where we have used n{t) = s{t) — h{t). Therefore the hkehhood ratio becomes 

p{s\h) p{s — h) 



A 



p(s|0) p{s) 
exp \—}^{s — h\s — h 



exp 



exp [— 



Ms) 



2{s\h)-{h\h)] + -{s\s) 



(4.30) 



exp 



[s\h)--{h\h) 



where {s\h) depends on the detector output and {h\h) is constant for a particular 
Sn{\f\) and h. Since the hkehhood ratio is a monotonicaUy increasing function of 
{s\h) we can threshold on {s\h) instead of the posterior probabilities. Our optimal 
receiver is involves the construction of {s\h) followed by a test 



[s\h) > x=K the signal is present, 
{s\h) < the signal is not present. 



(4.31) 



For a given /i(t), the inner product in equation (4.27), is a linear map from the 
infinite dimensional vector space of signals to M. Therefore the optimal receiver is 
a linear function of the input signal s{t). Both the output of a gravitational wave 
interferometer and inspiral signals that we are searching for are real functions of time, 
so 



r{f) = s{-f) 



and the inner product in equation (4.27) becomes 



(a I b) 



df 



(4.32) 
(4.33) 



(4.34) 



Snilfl) 

If we receive only noise, then the mean of {s\h) over an ensemble of detector 
outputs is 

{{s\h)) = {{n\h)) 



Sum 



(4.35) 
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since {n{t)) = 0. The variance of {s\h) in the absence of a signal is 



{{s\hr) 



oo poo 



CO J —CO 
oo POO 



■CO J —CO 
oo POO 



4 

{h\h) 



dfdf 
dfdf 
dfdf 



SnmSnilfl) 

, m)n*{f'))h*{f)h{f') 

Sum Sum) 
j5„(i/'i)5(/-/')/^*(/)Mr) 



(4.36) 



CO J —CO 



where we have used the definition of the one sided power spectral density from equa- 

(4.37) 



tion (4.23). In the presence of signal and noise, then the mean of (s|/i) is 

((n + /i|/i)) = ((n) + /i|/i) = {h\K). 

We can also show that the variance of {^s\K) in the presence of a signal is 

([(.|/^)-(/.|/.)f> = ([(n|/.)]^> = (/.!/.). (4.38) 

Therefore the quantity (/i|/i) is the variance of the output of the optimal receiver, 
(s|/i), and we denote it by 

= {h\K). (4.39) 
Now suppose that the signal we wish to recover has an unknown amplitude, A. 



The above discussion holds with h{f) Ah(t) and, from equation (4.30), the likeli- 
hood ratio becomes 



A = exp 



A{s\h) - -A\h\h) 



(4.40) 



which is again monotonic in {s\h), and so our previous choice of optimal statistic and 
decision rule hold. Now we are ready to consider the case of a gravitational wave 



inspiral signal of the form given in equation (4.1). The likelihood ratio now becomes 
a function of 9 

A\e) = p{e) exp !^D-\s\A{t) cos [20(t) -9])- ^D-\h\h) | . (4.41) 

Now consider the first inner product in the above exponential. Using cos(0 — 6) = 
cos 6 cos (p + sin 6 sin 0, we may write this as 

{s\A{t) cos [2(j){t) - 6]) = cose {s\A{t) cos [2(f){t)]) + sin^ {s\A{t) sin [20(t)]) 

= X cos 6^ + ?/ cos 6* (4.42) 

= \z\ cos(<l> - 9) 
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where 



X 

y 

\z\ 
tan<l> 



\z\ cos 

1^1 sin<l> = (s|A(t)sin(20(t))) , 

y_ 

X 



(4.43) 
(4.44) 
(4.45) 
(4.46) 



(The notation \z\ will become clear later in this chapter.) To calculate the likelihood 
ratio, A, we assume that the unknown phase is uniformly distributed between and 
2vr, 

p{d) = (4.47) 



and integrate A' over the angle 9 to obtain 



A 



27r 



M{e) do 



1 

2^ 



27r 



exp 



^" JQ 



D |2;|cos($ 



D- 



-{h\h) 



dd 



(4.48) 



where Jq is the modified Bessel function of the first kind of order zero. Once again, 
we note that the function /o(-D~^|2;|) is a monotonically increasing function of \z\ and 
so we can threshold on instead of A. Note that s appears in the expression for the 
likelihood through \z\ only. 

Recall from chapter [2] that we denoted the two orthogonal phases of the binary 



inspiral waveform by h^. and hg given by equations (|2.102|) and (|2.103|) 

/i 



hc{t) 
hs{t) 



2 



Mo 
ji_ 

Mo 



[7rGM/(t)]5cos [20(t) -20o] 
[7rGM/(t)]5sin [20(t) - 20o] 



(4.49) 
(4.50) 



and so for inspiral waveforms we can compute by 

z = ^{s\hcY + {s\hsY. 



(4.51) 



The threshold on |z| would be determined to achieve a given false alarm probability. 
We note that in the absence of signal |zp is the sum of squares of two independent 
Gaussian random variables of zero means and variance = (/^d^c) = (^sl^s)- ^ and 



65 



y are independent random variables since {hc\hs) = ^. It is therefore convenient to 
work with a normahzed signal-to-noise ratio defined by 



U|2 

•2 



/ = ^ (4.52) 



(t2 



which is distributed with two degrees of freedom for Gaussian detector noise. 

If a gravitational wave signal is present, then its location in time is defined by the 
end time parameter of the waveform. In chapter |2] we defined the end time of the 
chirp to be the time at which the frequency of the gravitational wave reached /ism, 
taken as the gravitational wave frequency of a particle in the innermost stable circular 
orbit of Schwarzschild spacetime. In the above discussion of the optimal receiver, we 
implicitly knew the location of the signal in the data to have = 0. Now suppose 
that the inspiral waveform ends at some unknown time te- We may write the signal 
we are searching for as h{t' — tg). Consider the Fourier transform of this signal 



e 



-2-/<'/i(t' _ Q dt' = e^^^'f'^ / e-^^'^^h{T) dr 



(4.53) 

where we have used t = t' — t^, dt = dr and t' = tg + t. The value of the inner 
product {s\hc) for a waveform that ends at time t^ is therefore 

{s\hM) = 2 ^/^''^^'^^^^fllf (4.54) 

and the signal-to-noise ratio for a chirp that ends at time t is 

p{t) = ^^{s\Kit)y + {s\K{t)y (4.55) 
a 

where the quantities {s\hc{t)) and {s\hs{t)) can be obtained by inverse Fourier trans- 



forms of the form in equation (4.54). 

Now the statistic p{t) derived from the likelihood is a function of a time parameter. 
For Neyman-Pearson optimal detection, we would integrate over all possible arrival 
times and threshold on this value. However, as well as making a statement about the 
presence or absence of a signal in the data we also want to measure the time that 
the signal occurs. To do this, we use the method of maximum likelihood[64\. The 
maximum likelihood estimator states that the most probable value for the location of 
the signal is the time at which the likelihood ratio is maximized. So to find a single 
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inspiral signal in a segment of interferometer data, we search for the maximum of 
p{t). If maxt [p{t)] > p* then we decide that we have detected a signal at the time 
of the maximum. When there is more than one inspiral in the data segment the 
maximization is not over all times. 

We have now completely specified the solution to the problem of finding a wave- 
form of unknown amplitude and phase at an unknown time in the data; our optimum 



receiver is the matched filter of equation equation (4.54). Below we develop the 
formalism to construct a digital implementation of the matched filter to search for 
gravitational wave signals in interferometer data. 

4.2 Conventions for Discrete Quantities 

The raw output of the interferometer is the error signal from the length sensing 
and control servo, LSC-AS_Q, as described in chapter [2j Although this signal is a 
dimensionless quantity, we say that it has units of "counts" and we denote it by v{t). 
The calibrated detector output is related to the raw detector output by the detector 
response function according to 

~s{f) = R{f;t)d{f) (4.56) 
where R{f;t) is the (complex) response function of the detector at time t and has 



units of strain/count (see section 2.2.3). In practice, the interferometer output is 
a discretely sampled quantity with sampling interval At, that is vj = v{tj) where 
tj = j At. The digital matched filter operates on a single data segment consisting of 
N consecutive samples of v{tj). The length of this data segment is T = NAt seconds. 
Henceforth, we let be a power of 2 and follow the convention that the subscript j 
refers to discretely sampled time domain quantities and the subscript k to discretely 
sampled frequency domain quantities. The frequency domain quantity v{fk) denotes 
the value of the continuous function v{f) at a particular frequency, labeled fk = 
k/{NAt). If the units of Vj are counts, then v{fk) has units of counts/Hz. We define 
the quantity Vk by Vk = v{fk)/At, which has units of counts. If k is negative, this 
corresponds to negative frequencies. 
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4.2.1 The Discrete Fourier Transform 

If v{tj) is sampled at intervals of At, then the sampling theorem [65] tells us that v{tj) 
is bandwidth limited to the frequency range —fny < / < Inj, where 

hy = ^ (4.57) 

is the Nyquist critical frequency. Any power in v{t) at frequencies above f^y will be 
aliased into the range — /nj < / < /Ny, corrupting the signal. To prevent this, signals 
of frequency higher than /nj in the interferometer output are removed using analog 
low-pass filters before the signal is digitized. Therefore v{tj) completely determines 
the signal v{t) in the band of interest. We may approximate the Fourier transform of 
this band limited signal v(tj) by 

N-1 N-1 

v{fk) Atv{tj)e-^^'f^'^ =^tYl ^je"'"''^/"^, (4-58) 

j=0 j=0 

where — (A^/2 + 1) < k < N/2 and the symbol means equal to under discretization. 
Notice that the approximation to the Fourier transform is periodic in k with period 
and so 

V-k = VN-k k = l,...,N -1. (4.59) 

Thus we let k vary from to — 1 where zero frequency (DC) corresponds to = 0, 
positive frequencies < / < /nj to values in the range < k < N/2 and negative 
frequencies — /nj < / < correspond to values in the range N/2 < k < N. The value 
k = N/2 approximates the value of the Fourier transform at both — /nj and f^y] both 
these values are equal due to the periodicity of the discrete transform defined by|66] 

N-1 
j=0 

We may estimate the discrete inverse Fourier transform in a similar way, using the 
relation 

A/ = Z^+i - ^ = 7^ - ]^ = ]^ (4-61) 

to obtain 

fc=0 
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4.2.2 Power Spectral Densities 



In equation (4.23), we defined the one sided power spectral density Sn{\f\) of n{t) to 
be 

WK(/)) = ^5„(|/|)5(/-/') (4.63) 

where angle brackets denote an average over different realizations of the noise. If n{t) 
has units of U then h{f) has units of (time) x U. The units 6{f — f) are (time), since 



5{f)df 



(4.64) 



is a dimensionless quantity and df has units (time)~^. Therefore we see that 5'„(|/|) 
has units of (time) x f/^. If we replace n{fk) with the discretely sampled quantities 

f^k = n{fk), we obtain 



N 
2At 



'S'n (|/fc|) ^kk' 



(4.65) 



where 6kk' is the dimensionless Kronecker (5-function, obtained by discretization of 
the continuous 5-function: 

6{f - /') ^ NAt6kk' (4.66) 



Equation (4.65) defines S'„ (|/fc|)in terms of the discrete frequency domain quantities. 
The definition in equation (4.65) is equivalent to 



Snilfkl) 



At 
N 

At 



k = 
k 



N 
2 ' 



(l^iV/2p) 

(d^feP + I^Af-feP)) otherwise 



(4.67) 



where the normalization is chosen so that the power spectral density satisfies the 
discrete form of Parseval's theorem 



N-l 
j=0 



N/2 



(4.68) 



k=0 



Parseval's theorem states that the total power in a signal is independent of whether 
it is calculated in the time domain or the frequency domain. 
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The value of Sn (|/fc|)for white Gaussian noise will be useful to us later, so we 
compute it here. If the noise n{t) is zero mean, white noise with variance then 



N-l N-1 
= ^ '^^27ri{jk-j'k')/N 

j=0 j'=0 

N-l N-l 
.Y^Y^^2.rUk-fk')/N^2^^^, 

j=0 j'=0 

N-l 

. ^ ^27rij(k-k')/N^2 
j=0 



Substituting this into equation (4.65), we obtain 

^ Sn (|/fc|) Skk' = N6kk'(^^ 



2At 



(4.69) 



(4.70) 



and so the power spectrum of white Gaussian noise is a constant with value 

SU\h\) = 2Ate. (4.71) 



4.3 Digital Matched Filtering 



The signal-to-noise ratio (4.52) requires us to compute the time series 

(/)C(/ 

Snilfl) 



x{t) = 2l d/e^-Z^^i^^lMZl 



and 



y{t) = 2l rf/e2-/*!i:^lMZ) 



(4.72) 



(4.73) 



Sum 

and the normalization constant a that measures that "amount of noise" in the detector 
(for a given inspiral waveform). From the definition of the inner product in equation 
(4.34) and the definition of o"^ in equation (4.39), we explicitly write 



oo 



oo 



(4.74) 



5.(1/1) J-oo Sj,i\f\) • 

The signal-to-noise ratio is normalized according to the convention of Cutler and 
Flanagan [27], so that in the case when the detector output is Gaussian noise, the 
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square of the signal-to- noise ratio averaged over an ensemble of detectors with different 
realizations of the noise is 

(P^> = ^(x^ + y^> = 2, (4.75) 



as seen from equation (4.36). 



4.3.1 Construction of the digital filter using stationary phase chirps 



In section 2.3.2 we derived the stationary phase approximation to the Fourier trans- 



form of the restricted post^-Newtonian binary inspiral waveform to be 

Un = 7^ (J^) ' (^) ' f-i f^V' e"<»«. (4.76) 



;iMpc)c2 V96M0y \ti^MqJ ' \ c3 
hsU)=^hcU), (4.77) 

where / is the gravitational wave frequency in Hz, M = mi + m2 is the total mass of 
the binary measured in solar masses, fi = mim2/M is the reduced mass and rj = fi/M. 
Note that hc,s{f) have units of 1/Hz and we have chosen the chirp to be at a canonical 
distance of r = 1 Mpc. The instrument strain per Hz h{f) is a linear superposition 
of hc,s{f) in the same way as h(t) is obtained from hc^s(t). The phase evolution to 
post^-Newtonian order is given by 

3 



^(/;M,r/)=27r/t,-2</)o-7r/4 + 



128r/ 

15 293 365 27145 3085 



5 /3715 55 , 3 2 



+ I \ ri H r?^ x"^ 

' 508 032 504 ' 72 ' y 

(4.78) 

where x = (nMfG/c^Y^^. The coalescence phase (po is the orbital phase, determined 
by the binary ephemeris, and the coalescence time tc is the time at which the bodies 
collide. The overall value coalescence phase (po is part of the unknown phase of 
the matched filter and we set 0o = 0, respectively. We set the coalescence time 



tc = 0, since it is accounted for by the Fourier transform in equations (4.72) and 



(4.73). The validity of the stationary phase approximation for inspiral templates is 
well established P^. 

Since the two chirp waveforms he and hg are orthogonal, the most efficient algo- 
rithm for constructing the time series p{t) uses a single complex inverse FFT rather 
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than computing it from x{t) and y{t) which requires two real inverse FFTs. We may 
further increase efficiency when using stationary phase chirps by sphtting the filter 
into a part that depends on the data and a part that depends only on the template 
parameters. In this section we describe the construction of a digital matched filter 



which uses these two tricks. Consider the discrete form of equation (4.72) 

7V-1 



_ o ^ \^ 27Tijk/N ^ifk)Kifk 



(4.79) 

At 2-Kijk/N ^kh*k 



where h^k = hc{fk)/ At. From equation (4.73) we obtain 



Recall that s(t) and h(t) are real signals. We may use the relations = /) and 



h{f) = h*{—f) to write the normalization constant, a^, defined in equation (4.74), as 

2 o 1 ^hc{fk)K{fk) 

k=0 
N-l 7 



At ^r~^ hckh*k 




Sum) S^ilfkl) S^lfkl) 

Since earth based gravitational wave detectors have no useful low frequency response, 
henceforth we set the DC {k = 0) term to zero. In addition to this, we assume that 
there is no power at the Nyquist frequency, as the low pass filter that band limits the 
interferometer data to frequencies below /nj falls off rapidly as the Nyquist frequency 
is approached. Therefore we may also set the k = N/2 term to zero. We assume this 
for all frequency domain quantities. 
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Now we may write the cosine phase of the filter given in equation (4.79) as 
At 



= 2— 



N-l 

E ■ 

k=N/2+l 
'n/2-1 



^2nijk/N _ 



N/2-1 



Sum) 



^2-Kijk/N ^kh, 



ck 



E 

k=l 



-2mjk/N ^k^ck 



k=l 



N/2-1 



Snilfkl)^ ^ ' 



Snilfkl) 

2mjk/N 



k=l 



Snilfkl) 



(4.82) 



where we have used the fact that fk = fN-k- Qj is defined to be 

Skkk 



N/2-1 



k=l 



SnilfklY 



(4.83) 



The sine phase of the filter given in equation (|4.80|) can similarly be written as 
At 



N 
At 



iV-l 

E ■ 

k=N/2+l 
'n/2-1 



^2TTijk/N ^kh*k 



Snilfkl) 
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_^ ^ ^2-Kijk/N 



SkK 
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E -2nijk/N ^l^sk , : 
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N/2-1 



2mjk/N ■^kii'sk 



Snilfkl) 

Skh 



k=l 



Snilfkl) 



k=l 



Snilfkl) 



(4.84) 



Using hs = ihc, equation (4.84) becomes 
At 



N 



■2nijk/N ^k^i^ck \ ^ ^2-nijk/N ^k\ '^j'T'ck 

Sn(IAI) ^ 



E 



k=l 



Snilfkl) 



At 



N/2-1 



N/2-1 



E -2TTijk/N ^k^ck \ ^ 2mjk/N ^kh*^k 

^ ^ Snilfkl) 



k=l 



k=l 



(4.85) 



2^^iQ*-Q,). 



Thus the outputs of the filter for the two phases are 



Xj 



(4.86) 
(4.87) 
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The quantity Zj is defined to be 



A ~ I* 

' N SA\fk\) 



At 



k=l 
N-l 



Y^^2.^,k|N~^ 



k=0 



where 



S„{|/fc|) U ^ ft. ^ 2 ' 



Zk = { (4.89) 
otherwise. 
We can now compute the square of the signal-to-noise ratio 

P%) = ^ = ^k/ (4.90) 

by a single complex inverse Fourier transform and threshold on > pi. Since we 
choose the template hc{f) to be at a canonical distance of 1 Mpc, the effective distance 
"D to a chirp detected with signal to noise ratio can be established as 

I?=-Mpc. (4.91) 
P 

Recall that a is a measure of the noise in the interferometer output; it is a measure 
of the sensitivity of the detector. Larger values of a correspond to a quieter detector 
(due to the 1/Sn{\fk\) term in the expression for a) and smaller values to a noisier 
detector. 

4.3.2 Details of Filter Implementation 

The calibrated detector output is related to the raw detector output by the detector 
response function according to 

~s{f) = R{f;t)d{f) (4.92) 

where R{f;t) is the (complex) response function of the detector at a specific time t, 
as described in chapter [2j In practice, we compute the uncalibrated power spectral 
density 5*^(1/^1) from the raw data and then the calibrated power spectral density, 



Sn{\fk\), in the denominator of (4.88) is 

Sn{\fk\) = \R{f;t)\'S.{\M). (4.93) 
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Further details of the computation of 5*^,(1/^1) are given in sections 4.6 and 4.7 



Typical values of the variance of v{t) for initial LIGO data are 10^; however, the 
response function R{f) has magnitude ~ 10~^^ at the most sensitive frequencies of the 
instrument. This means that Sn{\fk\) ~ 10"'*^ which is beyond the range of 4-byte 
floating point number^ We may store such values as 8-byte floating point numbers, 
but this is wasteful of memory since the extra precision of an 8-byte number is not 
needed. Therefore when we implement the digital filter, we multiply the response 
function R{f) by a scaling variable d which typically has values of d = 2^^ for initial 
LIGO data. This scales all frequency domain quantities to have approximately order 



unity. Therefore equation (4.88) becomes 

N/2-1 



k=l 



^2TTijk/N 



dRvk dh 



ck 



N 



d^\R\^S,{\M) 



and (4.81) becomes 



A N/2-1 

At 



d'^hckh*^k 



(4.94) 



(4.95) 



N ^ d^\R\^S.{\h\y 

Notice that we must also multiply the chirp by d so that all the factors of d cancel 
in the signal-to-noise ratio p{tj) and the normalization constant o"^. In fact this is 
convenient as it brings the value of the hc{fk) to order unity for chirps that would 
produce a signal-to-noise ratio of order unity. From equation (4.76) we obtain the 
dimensionless quantity 



d h, 



ck 



5fi 
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At 
2dGMQ 

(lMpc)c2 V96M0 
xexp[t^{fk;M,7])]Q{k-k 
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GMr, 







ttWo / V c^At 



(4.96) 



where the term Q {k — k^sco) ensure that the chirp is terminated at the frequency 
of the innermost stable circular orbit of Schwarzschild. The function \E'(/fc;M, r/) 
is the value of the post^-Newtonian phase evolution, which is given by (4.78), at 



^The smallest non-zero value that can be stored in an IEEE 754 floating point number is 
1.17549435 x IQ-^s. 
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the frequency fk- The quantity T{M,fi) in equation (4.96) is called the template 
dependent normalization constant and is given by 



T(M,/i) 



2dGMc: 







;iMpc)c2 J V96Mc 



5/i 



M 







n 2 



(4.97) 



Note that T{M, rf) depends on the masses of the template and as such must be 
recomputed once per template. If we substitute equation (4.96) into equation (4.95) 
we obtain 



a'=A 
where S is defined to be 



At 



^is 



At 



(4.98) 



(4.99) 



^-^^d^\R\^sA\fk\y 

S is referred to as the segment dependent normalization. It depends on the binary 
masses only through fcisco, so we compute and store the array 

N 



1 ^ k\^rn ^ 



(4.100) 



from the input power spectral density. We then select the correct value of S for a 
given mass pair by computing h^^o = /isco/A/. 
The signal-to-noise ratio squared is then 

16 fAtV 



\ N J 



r 



k=l 



d^\R\^S,i\h\) 



(4.101) 



where cr^ is now given by equation (4.98). Let us define qk by 

< k < 



Qk 



rf2|K|254|/fc|) 



(4.102) 







otherwise. 



and qj as the discrete complex inverse Fourier transform of qk- Then the signal-to- 
noise ratio squared is 



p\t, 



16 /AT 



(4.103) 



The computation of qk can further be split into the template independent computation 



of 



dvki^y 



d^\R\^S,{\fk\) 



(4.104) 
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and the computation of 

fk = exp [W{fk; M, 7])] e{k- fcisco) (4.105) 

where Fk is called the findchirp data segment and is called the findchirp template, 
so 

qk=\ ' ' (4.106) 

I otherwise. 

The goal of this separation is to reduce the computational cost of producing p(t) 
by computing the template and using it to filter several data segments. For a given 
power spectral density Sjj{\fk\) we only need to compute iS(A;isco) once. The findchirp 
code is designed to process several data segments, labeled i = 1, . . . , M, at a time. We 
can compute once for each data segment and then for each template we compute 
T and T^. This reduces the computational cost of filter generation. Furthermore, we 
can threshold against the quantity 

\q\l = ^lii^ (4.107) 

\ N J ^ 

thus saving a multiplication per sample point. The effective distance of an inspiral 



signal at time tj = jAt is given by equation (4.91) which becomes 



V . - (4.108) 



in this notation. 



4.3.3 Recording Triggers 



We call times when the optimal receiver tells us that a signal is present inspiral 
triggers and record the time of the trigger, the mass parameters of the template and 
the value of cr^ for the data segment. There are several complications that mean that 



simply thresholding on equation (4.107) is not what we do in practice, however. In 
section 4.4.2 we will show that an impulse in the data segment can cause the filter 
output event though no chirp is present, and hence cause a false alarms. Although 
such events are rare in Gaussian noise, they are quite common in real detector output, 
so we construct an addition test on the presence of absence of the signal, called the 
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veto|67], which is described in section 4.8 Furthermore, the inspiral signals that we 
are searching for are shorter than the length of a data segment, so we want to allow 
the possibility of generating multiple inspiral triggers in a single data segment. We 
also do not record all times for which Ig^-p > \q\l, as we would soon be flooded with 
triggers in real interferometer data. The algorithm that we use to select the times for 
which we generate inspiral triggers based on the output of the matched filter and the 



X veto is described in section 4.9 



4.4 Testing the filtering code 



4.4.1 Normalization 



Consider the case when the filter input is Gaussian noise, i.e. = rij. and set 
R{fk) = 1 and d = 1. Then the expectation value of the signal-to-noise ratio squared, 
(p2) is 



16 Z' ^ ^ ^isco ^isco 
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(4.109) 



where we have used the definition of a from equation (4.98) and the definition of the 



one-sided power spectral density from equation (4.65). 



The first test of the code is to check that the normalization of the filter agrees with 



equation (4.109) when the response function, i?, and dynamic range scaling, are 



both set to unity. In order to exclude issues related to power spectral estimation at 
this stage of testing we set the power spectral density to be the (constant) theoretical 
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value for white Gaussian noise given by 

Sni\fk\) = '2eAt, (4.110) 

where is the variance of the Gaussian noise. We generate five data segments 
containing white Gaussian noise of mean zero and variance q"^ = 64 at a sample 
rate of 16384 Hz. The length of each segment is 1920 seconds, so there are N = 
31457280 samples per segment. Table [l] shows the value of (p^) after averaging the 
output p^(tj) of the filtering code over all output samples. The values obtained are in 
good agreement with the theoretical expectation. Similar tests were performed with 
colored Gaussian noise, where the power spectrum is no longer a constant, and noise 
colored by a response function R{f); the average filter output was consistent with 
the expected value. Large and small values of the variance for the noise, were also 
used to test that the dynamic range scaling factor d was correctly implemented. In 



all cases the output of the filtering code was consistent with equation (4.109). 

We may also consider the distribution of the signal-to-noise squared in the presence 
of Gaussian noise. It is we can see from the definition of the filter x{t), given by 



equation (4.72), that it is a linear map from s{t) to x{t), and similarly for the filter 



in equation (4.73) that maps s{t) to y{t). If the input signal is s{t) is a Gaussian 
random variable, then the filter outputs x{t) and y{t) will be (uncorrelated) Gaussian 
random variables. Since the filter output p^{t) is the sum of the squares of these two 
Gaussian quantities, it will be distributed with two degrees of freedom. Recall 
that for a random variable, X, the cumulative density function is defined to be 

/X 

p{x)dx (4.111) 
-oo 

where p{x) is the probability density function. For a distribution with 2 degrees 
of freedom, this is 

Pix')= r ^dx. (4.112) 



Figure 17 shows the cumulative density function of p^(t) obtained from one of the data 
segments in Table [l] plotted against the theoretical value given in equation (4.112). 
Clearly the measured and theoretical values agree very well. 
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4.4.2 Impulse Time 

The second test is to examine tlie output of tlie filter in the presence of a delta 
function and a constant (white) power spectrum. The input to the matched filter is 

N-l 

h = Yl 5,7e-'"^''=/^ = e-2-'^V^. (4.113) 

k=0 



Substituting equation (4.113) into equation (4.101), we obtain 



p^t,) = hlit, - t,) + hlit, - t,) (4.114) 



which is the sum of the squares of the time reversed chirps. Figure [18] shows the 
output of the matched filter with a delta function input at t = 90 seconds. The 
length of the data segment is 256 seconds, the template has mi = m2 = IMq and the 
low frequency cut off of the template is 40 Hz. The length of this template is 43.7 
seconds. It can be seen that the filter output does indeed follow the form of equation 



(4.114). The impulse time is the time at which an impulse in the data would cause 



the filter output to peak at t^- We can see from equation (4.114) and figure 18 that 
for the filter we have implemented, the impulse time will he at t = to, since this is 
when the maximum of the filter occurs in the presence of an impulse. 

4.5 Wrap-around of the Fast Fourier Transform 

A simple experiment serves to demonstrate the effect of periodicity of the Fast Fourier 



Transform (FFT) in matched filtering. As with the example depicted in figure 18 



we generate an input data segment of length 256 seconds. Now we place the impulse 



at t = 250 seconds, however. Figure [20] shows the input and output of the filter for 
such a data segment. Notice that the output of the filter wraps around, so that the 
first 43.7 — 6 = 37.7 seconds of the filter output is non-zero. This is due to the Fast 
Fourier Transform treating the data as periodic: it identifies t = and t = 256. If 
the impulse was placed at t = 256, just before the end of the segment, then the first 
tc seconds of p(t) would be corrupted, where tc is the length of the chirp template. 
This demonstrates that data at the start of the segment is being correlated with data 
at the end of the segment due to the wrap-around of the FFT. This is obviously 
unphysical, so we consider the first tc seconds of the signal-to-noise ratio corrupted 
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and ensure that we do not consider this data when searching for inspiral triggers. We 



will return to this problem in section 4.7 when we consider the construction of the 



inverse power spectrum 1/Sn{\fk\) used in the filter. 
4.6 Power Spectral Estimation 

Interferometer data is not stationary over long periods of time, so we cannot simply 
compute a single value of 5'„(|/|) to be used in the matched filter for all time. We 
must use a power spectrum that gives the noise level at the time of the data segment 
that we are filtering. To do this we use Welch's method [68] to estimate the average 
power spectral density using data close in time to the segment we are filtering. 

A Welch power spectral density estimate is defined by an FFT length, overlap 
length and choice of window function. We require that the frequency resolution and 
length of the power spectrum are the same as those of the data Vk and template 
hck- If the data segment is of length points with a sampling interval of At, then 
the power spectrum must be of length N/2 + 1 points with a frequency resolution 
of A/ = l/(A^At). (It is possible to generate the average power spectral density at 
a different frequency resolution and then interpolate or decimate it to the correct 
frequency resolution, however.) For simplicity of implementation, the length of the 
data used to compute the power spectrum is the same as that used in the filter data 
segment. To construct the average power spectrum we take Aseg data segments of 
length A^ from near in time to the segment being filtered. Each segment overlaps its 
neighbors by Aoveriap sample points, so we need 

A^chunk = NX AT^eg - (AT^eg - 1) X A^^^erlap (4.115) 

input data points to compute the average power spectrum. The A'chunk input data 



points are called an analysis chunk. In section 4.3.1, we discussed filtering several 
data segments through each template in the filtering code; we will see later that the 
data segments used in the filtering code have the same length and overlap as those 
used to estimate the power spectrum. 

Recall that since we are computing a discrete Fourier transform of the input data, 
any power that is not at a sampled frequency in the power spectrum will bleed into 
adjacent bins. This is a particular problem for LIGO data where there are a lot 
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of spectral line features, caused by power line harmonics or mirror suspension wire 
resonances. These features contain a lot of power and, in general, their frequencies do 
not lie exactly at sampled frequencies. To prevent this power bleeding into adjacent 
bins, we apply a Hann window to the data before taking the Fourier transform. This 
is a standard technique and for further details we refer to the discussion in [65j. 

To construct an average power spectrum from the Nseg individual spectra that are 
computed, we average the A^seg values for each frequency bin. That is, the value of 
Sn{\fk\) at a frequency fk is the average of the Nseg values of the power spectra at fk- 
In the standard Welch computation of the power spectral estimate, the mean is used 
to average the values in each frequency bin. Consider using the mean to compute the 
average in the presence of a loud signal. If the data that contains the loud signal is 
used in the computation of the average spectrum, then Sn{\fk\) will contain power 
due to the signal. This will suppress the correlation of the signal and the template 
at those frequencies and cause the value of the signal-to-noise ratio to be lower than 
one would obtain if the average power spectrum is computed from noise alone. To 
avoid this problem we use the median to estimate the average power spectrum. This 
has two advantages: (i) computational simplicity, as we only need to compute one 
PSD and can use it for several data segments and (ii) insensitivity to outliers in the 
spectra, which means that excess power in one segment does not corrupt the spectra 
for neighboring segments. This is useful since the LIGO data is not truly stationary. 

For a Gaussian random variable, the median is a factor of log 2 larger than the 
mean. We must therefore divide the median power spectrum by log 2 to ensure that 
it has the same normalization as the mean power spectrum for Gaussian noise. This 
scaling has the unwanted effect of suppressing constant features in the spectrum, 
such as power lines and wire resonances, by a factor of log 2 compared with the mean 
spectrum. In practice we find that this does not have a significant effect on the output 
of the filtering code. For a low number of data segments, A^^geg, the log 2 correction 
factor is incorrect; the true value is between log 2 and 1 and we do not correct for this 



bias. Figure 21 shows the cumulative distribution of the filter output in the presence 
of Gaussian noise, where the average power spectrum is computed using the median 
method. The bias introduced by the for low N^^g does not have a significant effect on 
the filter output. 
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4.7 Computation of the inverse power spectrum 



We observed in section 4.5 that the FFT we use to compute the match filter treats 
the data as being periodic and that we had to ignore part of the filter output that 
was corrupted due to wraparound of the data. 



If we look at the correlation in equation (4.101), we can see that we are filtering 



the data against the inverse power spectrum as well as the chirp, that is our filter is 

Recall that the chirp has a duration that is typically much less than the length of the 
data segment, so the effect of wrap-around only corrupts a region that is the length 
of the chirp at the start of the data segment. Unfortunately, the length of the inverse 
power spectrum, as a time domain filter, is the same length as the data segment. 



Figure 19 shows the filter output when the input data is an impulse at t = 90 seconds 
and the power spectrum is computed from Gaussian noise using the median method. 
Notice that the filter output is non-zero at all times. No matter where the impulse 
is placed, the entire filter output would be corrupted by the inverse power spectrum. 
To prevent this, we truncate the square root of the inverse power spectrum to a 
length tinvspectrunc sccouds in the time domain. This means that the inverse power 
spectrum will have support (i.e non-zero values) for 2tinvspectrunc seconds in the time 
domain. Truncation of the inverse spectrum has the effect of smoothing out the 
high Q features (narrow line features, such as power line harmonics or resonances 
of the mirror suspension wires) and restricting the length of time that the filter is 
corrupted. The corrupted regions can then be ignored when searching for chirps in 
the filter output. 

The algorithm used to truncate the power spectrum is as follows: 

1. Compute the average power spectrum of the uncalibrated input data v{tj) using 
Welch's method as described in the previous section. 

2. Compute the square root of the inverse power spectrum. 



3. Set the Nyquist, {k = N/2) and DC {k = 0) components of this to zero. 
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4. Compute the inverse Fourier transform of \/S^^{\fk\) to obtain the time domain 
inverse PSD of length T = NAt seconds. 

5. Zero the square root of the inverse spectrum between the time tinvspectrunc/2 and 
(T — tinvspectrunc)/2 secouds. This sets the length of the square root of the inverse 
spectrum in the time domain to be tinvspectmnc seconds. 

6. Fourier transform the time domain quantity back to the frequency domain. 

7. Divide by the number of points to ensure that the inverse power spectrum is 
correctly normalized. 

8. Square this quantity to recover 5'^^(|/a:|)- 

9. Set the Nyquist and DC frequencies to zero. 

10. The (scaled) strain inverse power spectral density is then computed by 

,2c nf h = u Lm2 ^-''(I^^I)- (4.118) 

d^Sn{\fk\) \dxR{fk)\ 



The factor of 1/R{f) in equation (4.118) may add some additional length to 5^„(|/|) 
in the time domain (since R{f) is not white), but because R{f) is smooth with no 
sharp line features, this is insignificant. The length of the inverse power spectrum is a 
parameter that we may tune based on the nature of the data that we are filtering. In 
the analysis described in this thesis, we set the length of the inverse power spectrum 



in the time domain to 32 seconds. Figure 22 shows the filter output in the presence 



of an impulse for a truncated power spectrum. There is non-zero data before as well 
as after the impulse, so we must ignore data at the end of a segment as well as before. 

4.8 The veto 

Although the matched filter is very good at finding signals in the noise, transient 
events in the data will also cause high values of the signal-to-noise ratio for data 



containing impulses, as we saw in section 4.4.2 To distinguish a high signal-to-noise 



event due to a signal from one due to a transient, we use a time- frequency veto known 
as the veto. This was first proposed in [69] and is described in more detail in 
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In this section, we review the construction of the veto and in the next sections 
describe the implementation used in the filtering code. 

Let u and v be two orthonormal time series representing the two phases of a 
binary inspiral signal, hcitj) and hs(tj). We divide these waveforms into p frequency 
sub-intervals {ui} and {vi}, I = 1 . . .p with 



{ui\Um) = -Sim 

p 

{vi\Vm) = -Sim 
p 

{ui\Vm) = 

and u = X)f=i '^i and v = X)f=i ^i- 
We then obtain the 2p time series 

{xi} = {s\ui), 

{yi} = {s\vi), 

where s is the detector output. Notice that 

p 



X = = {h\u) 

1=1 
p 



1=1 

so that {x^ + y^)/cr^ is the signal to noise ratio squared p^. Now, let 

X 

Axi ^xi , 

p 

Ayi = yi — 
p 



4.119) 

4.120) 
4.121) 



4.122) 
4.123) 

4.124) 
4.125) 

4.126) 
4.127) 



and define 

= + (4-128) 

In the presence of Gaussian noise s — n this statistic is distributed with v — 2p — 2 
degrees of freedom. Furthermore, if a signal h = Au + Bv (with signal to noise squared 
of Psignai = A"^ + B"^) is present along with Gaussian noise s = h + n, then = is 
still distributed with v — 2p — 2 degrees of freedom. Small values of the veto 
mean that the signal-to-noise ratio has been accumulated in a manner consistent with 
an inspiral signal. We apply an additional threshold on for triggers that have a 
high signal-to-noise ratio. 



4.8.1 Implementation of the Digital Veto 

Recall that the templates h^. and hg are normalized such that 



and 



4 At ^ hck(i)h 



N 



k=0 



Sn{\fk\) P 



4At ^ 



W'^skim) 1 . 2 



Snilfkl) P 



fc(0'^sfc(m) 
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(4.129) 



We construct the p templates |/ic{o| |/is(/)|, where I = 1, . . . ,p, with 



(4.130) 
(4.131) 
(4.132) 



he = ^ hc{i), 
1=1 
p _ 

Z=l 



We construct the 2p filter outputs 



X 



and 



A:=0 



N/2 



SkK 

Sj]fk\) 



- A^sr^ .2.ijkiN ^^Kk{i) 



A:=0 



from which we can recover equations (4.79) and (4.80) by 



(4.133) 
(4.134) 



(4.135) 



(4.136) 



m 



1=1 



and 



i=\ 



(4.137) 



(4.138) 
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Consequently, the signal-to-noise ratio can be written as 

2 / p \ 2 



P 



Let 



and 



Ax 



'm 



X 



^ym = ym 



p 

y_i 
p 



and define the quantity 



P 



5^[(Ax,(o)'+(Ay,(,))^ 



1=1 



(4.139) 

(4.140) 
(4.141) 

(4.142) 



If at any time tj the signal-to-noise ratio exceeds the threshold p{tj) > p^, then we 
compute X^iij) for the data segment. We can then threshold on < X* to decide 



if the signal-to-noise event is consistent with a true inspiral signal. In section 4.8.2 



we discuss a modification to this threshold, if the template and signal are not exactly 
matched. 



4.8.2 Mismatched signal 

The waveform of an inspiral depends on the masses parameters M, r] of the two objects 
in the binary. If the template being used in the matched filter does not exactly match 
the true signal h'{M',ri') then the output of the matched filter will be smaller than 
if the template was correct. The mismatch can arise for any number of reasons; for 
example errors in the theoretical template mean that the post^-Newtonian waveform 
does not match the true inspiral signal (which becomes important at high masses, 
M > SMq) or errors in the calibration function R{f) may change the amplitude 
and/or phase of the signal in the data relative to the corresponding template. 

The loss of signal-to-noise ratio due to mismatch is accompanied by an increase in 
which requires a modification of the threshold. Suppose a signal. Aw, that is not 
exactly matched by m or i; is present in the data: s = Aw where w is the (normalized) 
and A is an amplitude. (Here we assume no noise.) With no loss of generality, orient 
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u and V such that {w\v) — 0. If w is nearly parallel to u, separated by some parameter 
difference Sx"' (which is small) then 

wc^u+ ^Sx'^ + - J x'^Sxf^ (4.143) 



so 



(.l„).^(.+li.."+l^.."..>) 



since (wjii) is a local maximum for w — u; thus (-^{u) — 0. Let us define the 
"mismatch" between w and u as 

ds^ = 1 - (^/;|xi) (4.145) 

and then 

ds' = ga,Sx-5x^ = -- [j^,\uj Sx-Sx^. (4.146) 

so 

Now we compute x^. We can ignore the v terms. Thus 

1=1 

= ^ E [(^1^')' - mu){h\ui)/p + {h\uf/p'] (4.148) 

i=i 

= ^ E(%')' - (%)'■ 

Using the Schwartz inequality {h\uiy < {h\h){ui\ui) we obtain 

{h\uif < {h\h) (4.149) 
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and so we may write 



< mh) - {h\uf] 

= [{w\w) - {w\uf] 
= A^[l-{l- ds^f] 
~ 2A^ds^. 



(4.150) 



Therefore, iih = Aw + n where w has a shght mismatch ds^ = 1 — {w\u), then has 
a non-central distribution with v = 2p — 2 degrees of freedom and a non-central 
parameter A = 2A^ds'^ and A^ = p^. 

A distribution with v degrees of freedom has a mean of v and a variance of 
2v in Gaussian noise; hence one often considers the quantity /v which would have 
a unit mean and a variance of two in the presence of Gaussian noise alone. A non 
central distribution with v degrees of freedom and non-central parameter A has a 
mean oi v + \ and a variance of 2(z/ + A) x [1 + A/(z/ + A)]. (The factor in brackets in 
the variance is always between 1 and 2, and is not really important for our purposes.) 
Thus, in this case, the quantity 'X^ jiv + A) has unit mean and variance of between 
two and four. 

For the reasons discussed above, a true signal will never be exactly matched by 
one of our template waveforms, so we wish to conservatively modify the threshold 
X* to allow for the case of a mismatched signal. While it is possible to construct 
constant confidence thresholds on the non-central x^ distribution for various signal 
events, a crude (but perhaps adequate) prescription is to threshold on the quantity 
X^/ (z/ + A), which is roughly equivalent to thresholding on x^/ {p + A?ds^^. Since real 
interferometer noise is not really Gaussian, it is not important to use the exact result 
for the non-central x^ distribution, though this could certainly be done. This choice 
of threshold is conservative as we would not reject signals more often than we expect 
if using the true threshold. In practice, we threshold on 



where p is the signal-to-noise ratio of the signal and 5^ is a parameter chosen to reflect 
the to be the largest amount mismatch that a true signal may have with the template 
waveforms. Since we do not know all the contributions to 5^ (in particular we do 
not accurately know the contribution from errors in calibration), we set 5^ by Monte 



x'<xl{v^p'n 



(4.151) 



Carlo techniques, which will be described in chapter [5] 
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4.9 Trigger selection algorithm 

The object of the search algorithm is to generate a hst of inspiral triggers. A trigger 
is a time at which there may be a binary inspiral signal in the data stream. The GPS 
time recorded in the trigger would correspond to the coalescence time of an inspiral 
signal, which corresponds to the time at which the signal-to-noise ratio squared is a 



maximum, as shown in figure 23 



We have seen in sections 4.5 and 4.7 that if the length of the chirp is tcMrp seconds 
and the length of the inverse power spectrum is tpso seconds then we must ignore 
^psD/2-|-tchirp seconds of data at the beginning of the data segment and tpsD/2 seconds 
of data at the end of the data segment due to wrap-around corruption. To simplify 
the data management, we ignore the first and last quarter of a segment; we test that 
^psd/2 -|- tchirp is less that one quarter of a segment and generate an error if it is 
not. The error informs the user that longer data segments must be used to avoid 
corruption of the filter output. 

The signal-to-noise ratio of a trigger must exceed the threshold pi and the 
statistic for the trigger must be less than the threshold value + P^*^^)- When 

generating triggers, we must consider the fact that the length of a data segment is 
greater than the length of a chirp, so may be multiple chirps in a single segment. 
We could simply examine the time series for sample points where > pi, however 
it is likely that for a true signal there will be many sample points above threshold 
corresponding to the same event. Similarly, if the data is noisy, we do not wish to 
generate a flood of events by considering every sample point above threshold a unique 
event. We address this by a trigger selection algorithm that we call maximization 
over a chirp. Figure |24] shows the algorithm for constructing the list of inspiral 



triggers. It can be seen from the algorithm in figure 24 that multiple triggers for the 
same template can be generated in one data segment. The coalescence times for the 
different triggers must be separated by at least the length of the template waveform. 

The list of inspiral triggers is the final output of the filtering code. For each trigger 
generated, we store the GPS time, mass parameters of the template waveform, the 
signal-to-noise ratio, the value of the veto, the effective distance V of the trigger 
in Mpc and the value of cr^. 
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The core of the inspiral analysis pipehnes that we construct is the generation of in- 
spiral triggers. Once we have generated the inspiral triggers from the matched filtering 
and code, we can test for coincidence between multiple interferometers, examine 
environmental data and auxiliary interferometer channels for associated artifacts, etc. 
Construction of an analysis pipeline is described in the next chapter. 
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Random Noise Generator Seed 




Var(p2(t)) 


7 


2.0118 


4.0312 


15 


2.0059 


4.0196 


19 


1.9965 


3.9911 


43 


1.9998 


4.0023 


69 


1.9936 


3.9846 



Table 1 : The mean and variance of the filter output p^{t) for five samples of white, Gaussian 
noise with a constant power spectrum of Sn (|/fe|) = 2<j^(5T. The observed values agree with the 
expected value for the mean and the variance, showing that the implementation of the matched 
filter is correctly normalized. 
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P^(t) 



Figure 17 : In the presence of Gaussian noise the expected filter output, is the sum of the 

squares of two Gaussian distributed quantities and so should be distributed with two degrees of 
freedom. This figure shows the cumulative distribution function (CDF) of the filtering code output 
and the expected analytic value. The filter input is white Gaussian noise of variance and a 
constant power spectral density of Sn {\fk\) = 2<i^5T. It can be seen that there is good agreement 
between the observed and expected values. 
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Figure 18 : The top panel shows the filter input which consists of an impulse at to = 90. The power 
spectrum is set to that of white uncorrelated noise. The bottom panel shows the output of the filter. 
The filter output is the sum of the squares of the time reverse chirps and the maximum of the filter 
output occurs at the time of the impulse. 
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Figure 19 : The top panel shows the filter input which consists of an impulse at to = 90. The 
power spectrum is computed from Gaussian noise of the same length of the input data using Welch's 
method. The bottom panel shows the output of the filter. Due to the fact that the duration of the 
inverse power spectrum l/5„ in the time domain is the same length as the data segment, the 
entire filter output is corrupted due to the wrap around of the FFT. 
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Figure 20 : The top panel shows the filter input which consists of an impulse at to = 250. The 

power spectrum is set to that of white uncorrelated noise. The bottom panel shows the output of 
the filter. The length of the chirp template is 43.7 seconds. Notice that the filter output is non-zero 
for the first 37.7 seconds of the output due to the wrap-around of the FFT. 
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P^(t) 



Figure 21 : In the presence of Gaussian noise the expected filter output, is the sum of the 

squares of two Gaussian distributed quantities and so should be distributed with two degrees of 
freedom. This figure shows the cumulative distribution function (CDF) of the filtering code output 
and the expected analytic value. The filter input is white Gaussian noise of length 256 seconds 
and the power spectrum Sni\Jk\) is computed from 15 segments of white Gaussian noise length 256 
seconds, overlapped by 128 seconds using Hann windowing and the median estimator. 
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Figure 22 : The top panel shows the input to the filtering code which is an impulse at t = 90 
seconds. The average power spectrum is computed from typical LIGO noise and then truncated to 
16 seconds in the time domain. The duration of non-zero filter output is also 16 seconds. 
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Input and Output of FindChirp Code 
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Figure 23 : Output time series from the filtering code for an inspiral chirp in the absence of noise. A 
(2.0, 2.0) m© inspiral chirp is generated using the post^-Newtonian time domain waveform generation 
and injected into the data. This is filtered using the post^-Newtonian stationary phase waveform. 
The signal to noise squared and time series are shown. The signal to noise squared is a maximum 
at the coalescence time of the template inspiral signal. This occurs slightly after the coalescence 
time of the injected signal. The difference in coalescence times is due to the different methods of 
generating the chirp signal. 
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Figure 24 : The algorithm used to generate inspiral triggers. For a given inspiral template we begin 
by calculating the length of the chirp and the filter output. For a data segment of length N, the 
first and last iV/4 points in the segment may be corrupted due to FFT wrap-around, so we ignore 
them. For the rest of the data segment, we step through the filter output looking for times when 
the signal-to- noise and threshold are satisfied (the main loop). If wc find a point that passes the 
threshold tests, we label it jmax and enter the maximization over chirp loop. This steps through the 
data looking for the any larger values of \qj\'^ within a chirp length (given by iVchirp) of the time 
imax- If a larger value of \qj\'^ is found, we reset jmax and keep looking for any larger values. If no 
larger value is found within a chirp length (or we reach the end of the uncorrupted data) we generate 
an inspiral trigger and save its time, mass, signal-to-noise ratio, value of the veto, number of x^ 
bin (p) and the value of for this data segment. 
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Chapter 5 

Detection Pipelines For Coalescing 
Binaries 



Chapter |4] described the algorithms that we use to generate inspiral triggers given an 
inspiral template and a single data segment. There is more to searching for gravi- 
tational waves from binary inspiral than trigger generation, however. To perform a 
search for a given class of sources in a large quantity of interferometer data we con- 



struct a detection pipeline. In section 5.1 we give an overview of the the components 
used in a pipeline and how they fit together. We then describe the building blocks of 
the pipeline in more detail. Section [5^ describes data quality cuts, which are used to 



discard data which is unsuitable for analysis. The application of trigger generation 



to the data is explained in section 5.3 The use of data from multiple interferometers 
is described in section 15. 4[ In section 15.51 we show how data from the interferometer 
that does not measure the gravitational wave signal can be used. Finally in section 



5.8 we describe the pipeline that has been constructed to search for binary neutron 



stars and binary black hole MACHOs in the S2 data. 
5.1 Construction of Inspiral Pipelines 



A detection pipeline is a sequence of operations that starts with the raw data from 
the interferometers and produces a list of candidate events. Figure 25 shows a simple 
pipeline to filter the data from a pair of interferometers labeled IFOl and IF02. We 
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only use data that comes from interferometers in stable operation. Running an in- 
terferometer is a complex process that requires human operators who are trained to 
lock the interferometers. Locking the interferometer is the process of bringing it from 
an uncontrolled state to a state where light is resonant in the interferometer. From 
a state in which the interferometer optics are freely swinging, the operators manu- 
ally align the optics of the interferometer using the interferometer control systems. 
They then direct the automated lock acquisition system |70j to bring the Fabry-Perot 
cavities into resonance with the beam splitter positioned so that the light at the anti- 
symmetric port is a minimum. The recycling cavity is then brought into resonance 
and the length sensing and control servo maintains the locked state by monitoring 
the motion of the optics and adjusting their position accordingly. In addition to the 
operator, members of the LIGO Scientific Collaboration trained in the operation of 
the interferometer are present at the observatory. The collaboration member on duty 
is known as the as scientific monitor or scimon. Once the interferometer is locked, 
the operator and scimon decide if the quality of the data being recorded is suitable 
to be flagged as science mode data. If the data is suitable it passes the first cut for 
gravitational wave analysis and the operators and scimons continuously monitor the 
data using various data monitoring tools. Lock is lost when the length sensing servo 
no longer has enough dynamic range to maintain resonance in the interferometer. 
This is typically caused by large seismic events which may be local to the observatory 
(e.g. liquid nitrogen tanks creaking as the expand in the sun) or of global origin (e.g. 
a large earthquake in China has caused loss of lock). Poor data quality in the inter- 
ferometer may also require a break in lock to remedy. Continuous locked operation 
has been maintained for up to 66.2 hours in the Hanford 4km interferometer. Seismic 
noise in the Livingston interferometer limited the longest lock to 6.93 hours during 
S2. 

It is possible that the operators or scimons may make mistakes in deciding that 
data should be flagged as science mode; they may forget to enable calibration lines, 
for example. There may also be noise sources or transient events in the data which 
make it unsuitable for analysis but which are not easily detectable in the control 
room while the data is being taken. To prevent such data from being used in an 
astrophysical analysis a list of data quality cuts is compiled. The manual selection of 
science mode data may be considered the first data quality cut. Additional tests of 
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data quality are described in section |5.2 



The chirp signals from compact binary inspiral depend on the masses and spins 



of the binary elements, as described in section 2.3 Searching for gravitational waves 
from BBHMACHOs requires signals which only depend on the masses mi and m2 of 
the binary. The inspiral signals from BBHMACHOs lie in some region of the template 
parameter space described by the component masses mi and m2. A single template 
is not sufficient to search for signals from the population, as the template for a given 
pair of mass parameters may not produce a high signal-to-noise ratio when used as a 
filter to detect a signal with different mass parameters. To search for signals from a 
region of parameter space we construct a bank of inspiral templates as described in 



section 5.3 The template bank is constructed to cover the parameter space in such 
a way that we do not discard any signals from our target population. 

We then use the bank of templates to filter the data for inspiral signals using the 
matched filter and veto discussed in chapter El This results in a list of inspiral 



triggers. For the pipeline shown in figure 25, a bank of templates is generated used 
to filter the data for inspiral signals for each interferometer. We describe how trigger 



generation is used in a pipeline in section in detail in section 5.3 



One of the most powerful methods of rejecting false alarms is coincidence between 
different interferometers. As described previously, there are three LIGO interferome- 
ters which are operated simultaneously during science runs. The HI and H2 detectors 
are co-located at the LIGO Hanford Observatory and the LI detector is located at 
the LIGO Livingston Observatory. A true gravitational wave should produce a signal 
in all operating detectors at the same time, up to the time delay for the arrival of the 
wavefront of the gravitational wave at the observatories. We therefore require time 
coincidence, which demands that inspiral triggers be present in all operating detec- 
tors simultaneously, with time offsets less than the light travel time between detectors 
plus the measurement error of detection. While coincidence is the most obvious use 
of multiple interferometers, other coincidence tests may also be used, e.g. demanding 
consistency of the waveform parameters between the triggers from two detectors, or 
consistency in the recovered amplitude of the signal relative to the sensitivity of the 
detectors. We describe these tests in section 15.41 

While the goal during data acquisition is to ensure that the data recorded is as 
stationary and Gaussian as possible, transient noise artifacts may still be present in 



103 



the data. For example, it has been seen that a person jumping up and down in the 
control room at the observatory will cause a burst of noise in the gravitational wave 
channel. There may also be occasional glitches in the interferometer control systems 
that cause a transient in the gravitational wave channel, despite the best efforts of 
the experimental team. To allow us to distinguish such events from true gravitational 
wave signals, we record several thousand data streams of auxiliary interferometer 
control channels and physical environment monitor (PEM) channels. Auxiliary chan- 
nels monitor the state of the servo loops that control the interferometer and include 
information about the pre-stabilized laser, the length sensing and control system and 
the input and output optics. PEM channels record the data from devices such as 
seismometers, magnetometers and microphones placed in and around the interferom- 
eter. These devices are designed to detect environmental sources that may couple to 
signals in the gravitational wave channel. This data can be used to construct vetoes 
of the inspiral triggers if a coupling can be identified between a noise source present in 
an auxiliary or PEM channel and inspiral triggers in the gravitational wave channel. 



We demonstrate this process with examples in section 5.5 



The final step in constructing a pipeline is to turn the various elements described 
above (data quality cuts, template bank generation, trigger generation, coincidence 
and vetoing) into code that can be executed in an automated way on a computing 
cluster. The execution of the code must ensure that all the input data has been 
analyzed and the components of the pipeline are executed in the correct sequence. 
We use a directed acyclic graph (DAG) to describe the work flow of a pipeline. For 
example, we may construct construct a DAG to execute the simple pipeline in figure 



25 on all the data from LI and HI recorded in S2. The DAG describing the work 



flow is submitted to a computing cluster via the Condor high throughput computing 
system [71] . The Condor DAG manager executes the pipeline described in the DAG 



that we generate. This process is described in more detail in section 5.8.3 



Implicit in the above discussion is that fact that there are many parameters that 
must be set at each stage of the pipeline. For example: What data quality cuts should 
we use? What signal-to-noise and thresholds should we use when generating the 
inspiral triggers? What coincidence tests should we apply and what should their 
parameters be? What auxiliary channels and PEM channels should be used as vetoes, 
and how do we apply these vetoes to inspiral triggers? Answering these questions is 
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the key to turning a pipeline into a full binary inspiral search; we call the process of 
selecting the parameters tuning the pipeline. In fact, pipeline tuning and construction 
of the pipeline are not entirely separate. After constructing a pipeline and initial 
tuning, we may decide to revisit the pipeline topology before performing additional 
tuning. 

When tuning the pipeline we may wish to minimize the false alarm rate, i.e. 
minimize the number of candidate events that are not due to inspiral signals. We may 
simultaneously wish to minimize the false dismissal rate to ensure that the pipeline 
does not discard triggers that are due to real signals. The false alarm rate can be 
studied by looking at candidate events in the playground. The false dismissal rate 
can be studied by injecting signals into the data, that is generating a known inspiral 
signal and adding it to the data before passing it through the pipeline. Injection of 



signals is described in section |5.7| and chapter [6j 

When tuning the pipeline for data that will be used to produce an upper limit, 
we must ensure that we do not introduce statistical bias. A bias in the upper limit 
could be introduced, for example, by selecting PEM channel events to veto particular 
inspiral triggers which were associated with the PEM events purely by chance. Clearly 
this could systematically eliminate true inspiral events artificially in a way that would 
not be simulated in the efficiency measurements described below. To do avoid such 
a possibility, we select 10% of all the data that we record as playground data. The 
data from GPS time [t, t + 600) is playground if 

t - 729273613 = mod(6370). (5.1) 

Playground data is selected algorithmically to provide a representative sample of the 
full data set. We are free to pursue whatever investigations we wish on the playground 
data. Although we do not use this data in the upper limit calculation, however we 
do not preclude the detection of a gravitational wave signal in this data. We describe 
the process of tuning the S2 binary black hole MACHO search in chapter [7| 

If we are using data from multiple interferometers we can measure the background 
rate of inspiral signals. We do this by introducing a time shift into the data from 
different detectors before passing it through the pipeline. If we assume that noise 
between the detectors is uncorrelated and the time shift is sufficiently large, as de- 



scribed in section 5.6, then any candidate events that survive the pipeline should be 
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due to noise alone and not astrophysical signals. By measuring the background rate, 
we can measure the false alarm rate of the pipeline which can be used for both tuning 
and the computation of the upper limit or detection confidence. 



5.2 Data Quality Cuts 

The theoretical matched filter is optimized for Gaussian data with a known, noise 
spectrum that is stationary over the time scale of the data analyzed. The filter there- 
fore requires stable, well-characterized interferometer performance. In practice, the 
interferometer performance is influenced by optical alignment, servo control settings, 
and environmental conditions. The list of science mode data provided by the opera- 
tors can contain times when the interferometer is not operating correctly. Unstable 
interferometer data can produce false triggers that may survive both the test and 
coincidence. Data quality cut algorithms evaluate the interferometer data over rel- 
atively long time intervals, using several different tests, or look for specific behavior 
in the interferometer to exclude science mode data that is unsuitable for analysis. 
To decide if we should exclude science mode data based on a particular data quality 
cut, we can examine the performance of the inspiral code on data which is flagged as 
suspect by the cut. 



5.2.1 Photodiode Saturation 



Figure 26 shows the signal-to-noise ratio and time series of the loudest candidate 
event that was produced by the LIGO SI inspiral search[T2]. Also shown is the 
filter output for a simulated inspiral with similar parameters that was injected into 
well behaved interferometer data. Notice that the time series of p(t), X^{t) and the 
raw data are very noisy around the time of the SI loudest candidate. In contrast, 
the time series for the simulated signal is very clean. On further investigation, it was 
determined that the photodiode that records the light at the anti-symmetric port had 
saturated at the time of the SI loudest event. The system that converts the light into 
an electronic signal for the length sensing and control servo had exceeded its dynamic 
range causing a noise transient in the data. The saturation is a symptom of poor 
interferometer performance. Photodiode saturations are caused by large bursts of 
noise in the gravitational wave channel which corrupt power spectral estimation and 
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matched filtering. A test was developed to monitor the gravitational wave channel 
for photodiode saturations and this test has become a data quahty cut for current 
and future searches. 



5.2.2 Calibration Lines 

A second example of a data quality cut is based on the presence of calibration lines, 
which were described in section |2.2.3 The calibration lines track the response of 



the instrument to mirror movement or a gravitational wave, which varies over time. 
Without the calibration lines, it is not possible to construct an accurate response 
function. Since the inspiral search needs correctly calibrated data, a simple data 
quality cut checks for the presence of the calibration lines in the data. If they are 
absent, the data is discarded. 



5.2.3 Data Quality Cuts Available in S2 

The full list of available data quality cuts and their meanings for S2 are show in 
Table [2j The table is divided into two sections, mandatory and discretionary data 
quality cuts. Mandatory cuts represent unrecoverable problems in data acquisition or 
calibration and so we must exclude these times from the list of science mode segments. 
Discretionary data cuts are optional for a particular search. For the inspiral search we 
decide whether or not to use a cut based on the performance of the trigger generation 
code in playground data when a particular data quality cut is active, as described in 
section [5.8.1 The times remaining after we apply data quality cuts are called science 
segments. 



5.3 Inspiral Trigger Generation 

Chapter |4] describes the algorithms that we use to determine if an inspiral from a 
binary of masses mi,m2 is present in a single data segment. The input to trigger 
generation is: 

1. The template, hck, drawn from a template bank. 

2. The data segment to be filtered, {vj} j G [0,A^], where vj is the raw (uncali- 
brated) interferometer output channel, LSC-AS_Q. A data segment is the unit of 
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analysis for the filter and is a subset of a science segment. 

3. An average power spectral density, Sv{\fk\), of the channel LSC-AS_Q, 

4. The instrumental response function, R{fk), which is required to calibrate the 
data and the power spectrum. 

In this section, we describe how these quantities are constructed and used to generate 
inspiral triggers. 

5.3.1 Template Banks 

The matched filtering described in chapter |4] has been used to detect the inspiral 
waveforms from binary neutron stars in the mass range 1 Mq < mi, m2 < 3Mq and 
binary black hole MACHOs in the mass range 0.2 Mq < 1712,1711 < IMq. Since 
each mass pair {mi, 7/12} in the space produces a different waveform, we construct a 
template bank, a discrete subset of the continuous family of waveforms that belong 
to the parameter space. The placement of templates in the bank is determined by 
the mismatch of the bank, M. The mismatch is the fractional loss in signal to noise 
that results when an inspiral signal, s, is not exactly correlated with a template in 



the bank, h. In terms of the inner product defined in equation 4.34 of section 4.1 it is 

M=l-— M^ . (5.2) 

Vih\h){s\s) 

If we consider the distribution of binaries to be uniform in space, then the fraction of 
events lost due to the mismatch of the template from a population is approximately 
M^, i.e the range is decreased by a factor of M. A mismatch of 3%, i.e. a 10% loss of 
event rate, is conventionally accepted as a reasonable goal for binary neutron stars. 
For binary black hole MACHOs, we reduce the minimal match of the template bank 
to 5%. This is due to the fact that the number of templates in the bank for a given 
interferometer noise curve scales as approximately fn^^l^, where rriyam = mi = m2 is 
the mass of the lowest mass equal mass template in the bank|72]. If we lowering the 
minimal match then the computational cost decreases as M~^. The loss in signal-to- 
noise ratio is balanced by the fact that we are searching for a population of binary 
black hole MACHOs in the galactic halo. This population is far from homogeneous 
with almost all signals expected to produce a signal-to-noise ratio far above threshold. 
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so this 5% loss in signal-to-noise ration will hardly constitute any loss in observed 
event rate. 

The scaling in the number of templates as a function of the lower mass of the 
bank parameter space is due to the fact that the number of templates required to 
achieve a given minimal match is a function of the number of cycles that the inspiral 
waveforms spends in the sensitive band of the interferometer. The more cycles the 
matched filter correlates against, the greater its discriminating power and so the loss 
in signal-to-noise ratio for a mismatched template increases. A pair of inspiralling 
1.4 M0 neutron stars have 347 cycles in the S2 sensitive band (between 100 Hz and 
2048 Hz), compared to 1960 cycles for a pair of 0.5 Mq binary black hole MACHOs; 
a pair of 0.1 Mq binary black hole MACHOs have nearly 28 500 cycles. 

Since the number of templates is a function of the number of cycles of an inspiral 
in the sensitive band of the interferometer, it also depends on the shape of the power 
spectral density of the noise curve. In fact location of the templates in the bank is also 
a function of the PSD, as described in [72]. The power spectrum of the instrument 
changes over time and so we must also change the template bank. We accomplish 
this by using the power spectral density for an analysis chunk to generate a template 
bank that is unique to that chunk. The PSD is calibrated and the template bank gen- 



erated. Figure 29 shows binary neutron star and binary black hole MACHO template 
banks generated for a typical stretch of S2 data. The smallest and largest number of 
templates in a bank during S2 was 589 and 857 binary neutron star templates and 
11588 and 17335 binary black hole MACHO templates. 

5.3.2 Data Management 

Corruption due to the wrap-around of the matched filter means that not all the 
time in a data segment can be searched for triggers. As described in chapter |4| we 
simplify the process of selecting uncorrupted data by only searching for triggers in 
the signal-to-noise ratio, p^itj), when 

N 3N , , 

4<J<^, (5.3) 

where is the number of sample points in the data segment, and by demanding that 
the amount of data corrupted is less than A^/4 sample points. To ensure that all data 
is analyzed we must overlap each data segment by {N/2)At seconds. To compute 
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an average power spectrum we require a sample of data near the data segment being 
filtered so that a good estimate of the noise can be obtained. We combine these two 
requirements by bundling several overlapping data segments together in an analysis 
chunk. The length of an analysis chunk is bounded above by the amount of memory 
available in the computers performing the filtering and bounded below by requiring a 
sufficiently large number of segments in the computation of the average power spec- 
trum. In the S2 pipeline, we construct analysis chunks of length 2048 seconds from 15 
overlapped data segments of length 256 seconds. The data segments are overlapped 
by 128 seconds, with the first and last 64 seconds of data segment ignored when 
searching for triggers. The analysis chunks themselves are overlapped by 128 seconds 
so that only the first and last 64 seconds of a science segment are not searched for 
inspiral triggers. 



Figure 28 shows how analysis chunks are constructed from the science segments 
in S2. The first analysis chunk is aligned with the start of the science segment. 
Subsequent chunks overlap the previous one by 128 seconds. At the end of a science 
segment there is generally not enough data to fit an entire analysis chunk but we 
cannot make the chunk shorter, as we need all 15 data segments to compute the 
average power spectrum. To solve this problem, we align the end of the last chunk 
with the end of the science segment and ignore any inspiral triggers generated for 
times that overlap the previous chunk. 

The interferometer records data at 16 384 Hz, and we down sample the analysis 
chunk after reading it from disk to decrease the computational resources required by 
the filtering code. We choose the new sample rate so that the loss in signal-to-noise 
ratio due to the discrete time steps. At, is less than that due to the discrete choices of 
the template mass parameters. It can be shown that for the initial LIGO noise curve, 
this is true if the sample rate of the filtered data is greater than ~ 2600 Hz [72]. For 
simplicity, we chose sample rates that are powers of two and so we resample the data 
to 4096 Hz. We do this by applying a finite impulse response (FIR) low pass filter to 
remove power above the Nyquist frequency of the desired sample rate, 2048 Hz. The 
low passed data is then decimated to the desired sample rate. Although the maximum 
frequency of most of the BBHMACHO inspiral signals is greater than 2048 Hz, the 
loss of signal-to-noise ratio above this frequency is negligible as most of the the signal- 
to-noise is accumulated at frequencies lower than 2048 Hz, as described in chapter 
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|2} Figure |27| shows the loss in signal-to-noise ratio of a BBHMACHO inspiral due to 
resampling. The inspiral waveform of a pair of 0.2 Mq black holes an effective distance 
of 25 kpc is generated at a sample rate of 16 384 Hz. The maximum frequency of 
this waveform is 10112 Hz. The waveform is injected into raw (un-resampled) data 
with a typical S2 noise curve which is then filtered at 16 384 Hz and 4096 Hz. The 
maximum of the signal-to-noise ratio is p = 73.67 for the raw data and p = 73.58 
for the resampled data, giving a difference in signal-to-noise ratio of 0.1%. This loss 
combines the effects due to the discreteness of the resampled data and the signal 
present above the Nyquist and is much less than the 5% loss in signal-to-noise ratio 
caused by the discrete nature of the template bank. 

The interferometer data contains a large amount of power of seismic origin at 
low frequencies. This power is several orders of magnitude higher than the noise in 
the sensitive frequency band of the interferometer. This power may bleed across the 
frequency band when the data is Fourier transformed into the frequency domain and 
dominate the true noise in the sensitive band of the interferometer, In order to prevent 
this, we apply a Butterworth infinite impulse response (HR) high pass filter to the 
resampled data. The high pass frequency and filter order are chosen so that power 
in the sensitive band of the interferometer is not attenuated. The data is filtered 
forwards and backwards through the filter to remove the dispersion of the filter. The 
data used to compute the power spectral estimate is windowed using a Hann window 
to prevent power from line features in the spectrum (e.g. power line harmonics or 
suspension wire resonances) from bleeding into adjacent frequency bands. We also 
apply a low frequency cutoff in the frequency domain at a slightly higher frequency 
than the time domain filter. The matched filter correlation is not computed below 
this cutoff, so frequencies below it do not contribute to the signal-to-noise ratio. 

A windowed copy of each of the 15 data segment is used to construct the average 



power spectral density used in the matched filter, as described in section 4.6 Note 
that the data used in the matched filter is not windowed; the windowed data is 
discarded once it has been used to compute the power spectra. The mean values over 
the analysis chunk of the calibration parameters a and f3 are used to construct the 



response R{fk), as described in equation (2.62) of section 2.2.3 The same response 
function is used to calibrate the power spectral density and all data segments in the 
chunk. 
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A disadvantage to the above method of processing the input data is that we require 
a science segment to be at least 2048 seconds long. Any shorter science segments are 
ignored, as there is not enough data to generate a power spectral density. We note here 
that this lead to a significant amount of data being discarded in the S2 analysis. As 
we will describe below, the S2 pipeline requires the LI interferometer to be operating 
in order to analyze the data. LI is the least stable of the interferometers as it is 
very sensitive to seismic noise. High seismic noise can saturate the length sensing 
and control servo and cause the loss of lock, which terminates a science segment. 
Modification of the data management (e.g. construction of analysis chunks and power 
spectral estimation) to allow us to use science segments shorter than 2048 seconds 
requires significant changes of the implementation of the inspiral search code. It was 
not possible to implement and test these changes within the time allowed to perform 
the S2 analysis. Fortunately, it is expected that after the installation of addition 
seismic isolation at the Livingston observatory, scheduled for completion in late 2004, 
the lengths of science segments will be significantly increased and the number of short 
segments discarded will decrease. Unfortunately data from the third science run, S3, 
which was completed before the seismic upgrade, exhibits the problem of short science 
segments, so a redesign of the filtering code may still be required. 

5.3.3 Trigger Generation Pctrameters 

The process of template bank construction and inspiral trigger generation relies on 
several parameters that can be tuned to minimize the false dismissal or false alarm 
rate. We have touched on some of these already; in this section we enumerate all of 
the tunable parameters for bank and inspiral trigger generation. 

Both template bank and inspiral trigger generation require a calibrated power 
spectral density. In the inspiral trigger generation described above, construction of 
the PSD is coupled to the length of the data segments and analysis chunks. The 
following data conditioning parameters are used to construct the data segments, and 
so determine the characteristics of the PSD: 

• Number of sample points in a data segment, N. This determines the length of 
the segments correlated in the matched filter and the length of the segments 
used in the PSD estimate. The number of points that subsequent data segments 
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are overlapped by is then Ai'overiap = N/2. 

• Number of data segments in a chunk A'^segments- This sets the number of data 
segments used in the PSD estimate and so the number of segments in an analysis 
chunk. 

• Sample Rate, 1/At. The sample rates used in LIGO data analysis are integer 
powers of two Hertz. 

• Number of non-zero points in the square root of inverse power spectrum in the 
time domain, A^invspectrunc- This parameter was described in detail in section 

We set At A/'invspoctrunc = 16 sccouds for the S2 analysis based on the length of 
wraparound corruption allowed. A systematic study of the value of A^invspcctmnc has 
not been carried out for the SI or S2 data, but is planned for future analysis. 

Once the above parameters have be specified, it follows that the number of sample 
points in an analysis chunk is 

A^chunk (A^scgmcnts l)A^stride + A^ (5.4) 

where 

A'stride = N — A^ovcrlap- (5.5) 

The length of the analysis chunk, in seconds, is therefore 

^chunk = A'chunk- (5.6) 

The choice of these parameters is governed by the class of waveforms searched for and 
the low frequency sensitivity of the interferometer. The longest inspiral waveform in 
the template bank (which will be the smallest mass template), is determined by the 
lowest sensitive frequency of the interferometer. The sum of the length of the longest 
template and the length of the inverse power spectrum must be shorter than the 
duration of the signal-to-noise output, p(t), that we ignore due to corruption. We 
therefore require 

^A'^ ^ 2A^iiivspectrunc ~l~ A^longest (5' 7) 

where A^iongest is the number of points in the longest chirp. 
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Other data conditioning parameters control the cut-offs of the time and frequency 
domain filters applied to the data. These are: 

• The high pass filter cutoff, /hp, and high pass filter order Ohp. These parameters 
determine the shape of the IIR Butterworth high pass filter applied to the 
analysis chunks before PSD estimation and inspiral trigger generation. 

• The frequency domain low frequency cut off, /low This parameter allows us to 
exclude frequencies below a certain value from the correlations in the matched 
filter and veto. A non-zero value of /low sets the value of the data in all 
frequency bins k < /cmin = /low/ (NAt) to zero which excludes this data from 
the correlation. Note that /low > /hp to prevent data used in the correlation 
being attenuated by the high pass filter. 

During investigation of inspiral triggers in the S2 playground data, it was discovered 
that many of the LI inspiral triggers appeared to be the result of non-stationary noise 
with frequency content around 70 Hz. An important auxiliary channel, LI : LSC-POB_I, 
proportional to the residual length of the power recycling cavity, was found to have 
highly variable noise at 70 Hz. There are understandable physical reasons for this, 
namely the power recychng servo loop (for which L1:LSC-P0B_I is the error signal) 
has a known instability around 70 Hz, which often results in the appearance of glitches 
in the detector output channel at around 70 Hz. As a consequence, it was decided 
that to reduce sensitivity to these ghtches the high pass cut off should be set to 
/hp = 100 Hz with order Ohp = 8 and the low-frequency cutoff set to /low = 100 Hz. 
This subsequently reduced the number of inspiral triggers (presumably created by this 
problem); an inspection of artificial signals injected into the interferometer revealed a 
very small loss of efficiency for binary neutron star inspiral and BBHMACHO signal 
detection resulting from the increase in the low frequency cutoff. 

After we have defined the data conditioning parameters and a parameter space 
for the search, the only remaining free parameter for template bank generation is: 

• The template bank mismatch, M e [0, 1). Given a value of mismatch, M, for a 
real signal that lies in the bank parameter space the fractional loss in signal-to- 
noise ration should be no larger than M when filtering with the signal against 
its exact waveform compared to filtering the signal against the template in the 
bank with yields the highest signal-to-noise ratio. 
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As described above, the value chosen for the S2 search is Mbbhmacho = 5%. Note 
that we do not truncate the power spectrum used for generating the template bank, 
as there is no issue with wrap-around in the bank generation algorithm. 

For a given template we use matched filtering to construct the signal-to-noise 
ratio, p, and search for times when this exceeds a threshold, p > p*. If the threshold 
p^, is exceeded, we construct the template based veto, x^; with p bins. Small values 
of indicate that the signal-to-noise was accumulated in a manner consistent with 
an inspiral signal. If the value of the veto is below a threshold, < + ^'^p'^)) 
then an inspiral trigger is recorded at the maximum value of {p|x^ < + ^"^P^)}- 
The parameters used available in trigger generation are: 



The signal to noise threshold, p*. 
The threshold, x^- 

The number of bins used in the veto, p. 



The parameter 5"^ used to account for this mismatch of a signal and template 
in the x^ veto, as described in section 4.8.2 



Tuning of these trigger generation parameters is particular to the class of search 
used. A detailed discussion of this tuning is for binary black hole MACHOs is given 
in chapter [7j 



5.4 Trigger Coincidence 

Coincidence is a powerful tool for reducing the number of false event candidates 
surviving a pipeline. The simplest test of coincidence is time coincidence of triggers 
between two or more interferometers. For a trigger to be considered coincident in 
two interferometers, we demand that it is observed in both interferometers within a 
temporal coincidence window 5t. The coincidence window must allow for the error 
in measurement of the time of the trigger. It must also allow for the difference 
in gravitational time of arrival if the interferometers are not located at the same 
observatory. The time difference between gravitational wave arrival time varies from 
seconds, if the gravitational wave is propagating perpendicular to a line joining the 
detectors, to 10 ms if the gravitational wave is propagating parallel to a line joining the 
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detectors. The maximum time difference comes from the time it takes a gravitational 
wave to propagate between the observatories, given by t = ^, where d = 3002 km is 
the distance between the observatories and c is the speed of hght (the propagation 
speed of gravitational waves). Monte Carlo analysis with simulated signals suggests 
that we cannot measure the time of the trigger to an accuracy of less than 1 ms. 
The time coincidence window is therefore St = 1 ms if the interferometers are located 
at the same observatory. For coincidences between LHO and LLO triggers, we set 
6t = 3000 km/c + 1ms = 11 ms. Any triggers that fail this test are discarded. 

If a signal found in temporal coincidence is generated by real inspiral, it should 
have the same waveform in both interferometers, up to issues of the different detector 
antenna patterns yielding different combinations of /i+ and h^. This suggests that 
we could apply a waveform parameter test to triggers candidate that survive the 
time coincidence test. We cannot exactly extract the parameters of the waveform, 
however, since we filter the interferometer data with a template bank which may not 
contain the true waveform. In addition to this, the template banks will, in general, 
differ between detectors and detector noise may cause error in the measurement of 
the signal parameters, even if the template banks are the same. To account for these 
sources of error, we can apply waveform parameter coincidence by requiring that the 
two mass mi and 1712, of the templates are identical to within an error of 6m. 

We now consider an amplitude cut on the signals. The Livingston and Hanford 
detectors are not co-aligned. There is a slight misalignment of the detectors due to 
the curvature of the earth and so the antenna patterns of the detectors differ. This 
causes the measured amplitude and phase of a gravitational wave to differ between 
the sites. In the extreme case, it is possible, for example, for a binary to be completely 
undetectable by the LI detector, but still detectable by the HI and H2 detectors. For 
a given inspiral trigger, we measure the effective distance of the binary system. This 
is the distance at which an optimally oriented binary would produce the observed 
signal-to-noise ratio in a particular instrument — it is not the true distance of the 
binary. Since the detectors have different antenna patterns they will report different 



effective distances for the same gravitational wave. Figure 30 shows the ratio of 
effective distances between the two LIGO observatories for the population of binary 
neutron stars considered in the S2 analysis. The significant variation of the ratio 
of the effective distances precludes using a naive test for amplitude coincidence. It 
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is possible to obtain information about sky position from time delay between sites 
to construct a more complicated amplitude cut, but this has not be used in the S2 
analysis. 

In the case of triggers from the HI and H2 interferometers that are coincident in 
time and mass, we can apply an amplitude cut that tests that the effective distances 
of the triggers are coincident. In this test we must allow for the relative sensitivity of 
the detectors while allowing for error in the distance measurement, as determined by 
Monte Carlo simulations. The amplitude cut for triggers from HI and H2 is given by 

' ' + (5.8) 



where Di {T>2) is the effective distance of the trigger in the first (second) detector and 
p2 is the signal-to-noise ratio of the trigger in the second detector, e and n are tunable 
parameters. In order to disable the amplitude cut when comparing triggers from LLO 
and LHO, we set k = 1000. When testing for triple coincident triggers we accept 
triggers that are coincident in the LI and HI detectors that are not present in the H2 
detector if the effective distance of the trigger is further than the maximum distance 
of H2 at the signal-to-noise ratio threshold at the time of the candidate trigger. Figure 



31 summarizes the algorithm for the time, mass and distance coincidence tests used 
in S2. 

We therefore have the following coincidence parameters that must be tuned for 
the pipeline: 

• The time coincidence window, 5t, which is set to 1 ms for LHO-LHO coincidence 
and 11 ms for LHO-LLO coincidence. 

• The mass coincidence window, 5m. 

• The error on the measured effective distance due to the instrumental noise, e, 
in the amplitude test. 

• The systematic error in measured effective distance, k, in the amplitude test. 

If coincident triggers are found in HI and H2, we can get an improved estimate 
of the amplitude of the signal arriving at the Hanford site by coherently combining 
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the filter outputs from tlie two gravitational wave channels, 



12 



/ \zhi + zh2\ ^^ 
PH = W 9 , o , (5.9) 



where z is the matched filter output given by equation (4.88). In this combination, 
the more sensitive interferometer receives more weight in the combined signal-to- 
noise ratio. If a trigger is found in only one of the Hanford interferometers, then pn 
is simply taken to be the p from that interferometer. 

Finally, we cluster the coincident triggers over a 4 second time interval. Clustering 
is performed so that a noise transient that may cause several templates to trigger 
within a small window is only counted as a single event in the data sample. 

5.5 Auxiliary and Environmental Channel Vetoes 

In addition to data quality cuts, another method to exclude false alarms is to look for 
signatures in environmental monitoring channels and auxiliary interferometer chan- 
nels which would indicate an external disturbance or instrumental glitches. This al- 
lows us to veto any triggers recorded at that time. Auxiliary interferometer channels 
(which monitor the light in the interferometer at points other than the antisymmetric 
port — where a gravitational wave would be most evident) are examined, with the aim 
being to look for correlations between glitches found in the readouts of these channels 
and inspiral event triggers found in the playground data. By doing so, we are capable 
of identifying instrumental artifacts that directly affect the light that is measured in 
the gravitational wave channel, so these vetoes are potentially very powerful. Figure 



32] demonstrates the the use of auxiliary channels to identify the source of an inspiral 
trigger in the gravity wave channel. 

When choosing vetoes, we must consider the possibility that a gravitational wave 
itself could produce the observed glitches in the auxiliary channel due to some physical 
or electronic coupling. This possibility was tested by means of hardware injections, 
in which a simulated inspiral signal is injected into the interferometer by physically 
moving one of the end mirrors of the interferometer. Hardware injections allow us to 
establish a limit on the effect that a true signal would have on the auxiliary channels. 
Only those channels that were unaffected by the hardware injections were considered 
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"safe" for use as potential veto channels. The process of testing veto safety with 
hardware injections is described in more detail in chapter |6| 

We used a computer program, alitchMonf73\ . to examine the data and identify 
large amplitude transient signals in auxiliary channels. Numerous channels, with var- 
ious filters and threshold settings, were examined and which produced a list of times 
when the glitches occurred. The glitch event list was compared with times generated 
by triggers from the inspiral search (Note that these studies were all conducted on 
the playground data.) A time window around a glitch was defined, and any inspiral 
event within this window was rejected. One can associate the veto with inspiral event 
candidates and evaluate a veto efficiency (percentage of inspiral events eliminated), 
use percentage (percentage of veto triggers which veto at least one inspiral event), 
and dead-time (percentage of science-data time eliminated by the veto). A "good" 
veto will have a large veto efficiency and use percentage with a small dead time sug- 
gesting that it is well correlated with events in the gravitational wave channel that 
produce inspiral triggers. Followup studies are performed on such candidate vetoes to 
determine the physical coupling between the auxiliary channel and the gravitational 
wave channel. If a sound coupling mechanism is found, then the veto will be used. 

Tuning of the vetoes for binary black hole MACHOs is described in chapter [7j 

5.6 Background Estimation 

Since we restrict the S2 analysis to coincident data and require that at least two 
of the interferometers must be located at different observatories, we may measure a 
background rate for our analysis. We estimate the background rate by introducing 
an artificial time offset, or lag, T to the triggers coming from the Livingston detec- 
tor relative to the Hanford detector. We call this "sliding the triggers by T." After 
generating triggers for each interferometer, we slide the triggers from the LHO inter- 
ferometers relative to the LLO interferometer and look for coincidences between the 
offset and zero lag triggers. The triggers which emerge from the end of the pipeline 
are then considered a single trial representative of an output from a search if no 
signals are present in the data. By choosing a lag of more than 20 ms, we ensure 
that a true gravitational wave will not be coincident in the time-shifted data streams. 
In fact, we use lags longer than this to avoid correlation issues; the minimum lag is 
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17 seconds. Note that we do not time-shift the two Hanford detectors relative to one 
another since there may be real correlations due to environmental disturbances. If 
the times of background triggers are not correlated in the two interferometers then 
the background rate can be measured; we assume that there is no such correlation 
between LHO and LLO triggers. 

5.7 Detection Efficiency 

In absence of detection, we will construct an upper limit on event rate. To do this 
we need to measure the detection efficiency, e, of the analysis pipeline to our pop- 
ulation. This is the fraction of true signals from a population that would produce 
triggers at the end of the pipeline. A Monte Carlo method is used to measure this 
efficiency. We simulate a population of binary neutron stars and inject signals from 
that population into the data from all three LIGO interferometers. The injection is 
performed in software by generating an inspiral waveform and adding it to interfer- 
ometer data immediately after the raw data is read from disk. We inject the actual 
waveform that would be detected in a given interferometer accounting for both the 
masses, orientation, polarization, sky position and distance of the binary, the antenna 
pattern and calibration of the interferometer into which this signal is injected. The 
effectiveness of software injections for measuring the response of the instrument to 
an inspiral signal is validated against hardware injections where an inspiral signal 
is added to the interferometer control servo during operation to produce the same 
output signal as a real gravitational wave. This validation is described in chapter |6| 
The data with injections is run through the full analysis pipeline to produce a list of 
inspiral triggers. We may combine the signal-to-noise rations from coincident triggers 
from several interferometers into a single coherent signal-to-noise ratio. 



where the form of / is chosen based on studies of the playground and background 
triggers. We can then construct a final threshold, p^,, on triggers that survive the 
pipeline. The detection efficiency, e{p), is the ratio of the number of signals with 
p > to the number of injected signals. 




(5.10) 
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5.8 The S2 Data Analysis Pipeline 

In this section we describe the pipehne constructed to search the S2 data for gravita- 
tional waves from inspiralhng binary neutron stars and binary black hole MACHOs. 
The data quality cuts used are common to both searches and are described in section 



5.8.1, The detection of a gravitational-wave inspiral signal in the S2 data would (at 
the least) require triggers in both LI and one or more of the Hanford instruments with 
consistent arrival times (separated by less than the light travel time between the de- 
tectors) and waveform parameters. During the S2 run, the three LIGO detectors had 



substantially different sensitivities, as can be seen from figure 33, The sensitivity of 
the LI detector was greater than those of the Hanford detectors throughout the run. 
Since the orientations of the LIGO interferometers are similar, we expect that signals 
of astrophysical origin detected in the Hanford interferometers will most often be also 
detectable in the LI interferometer. We use this and the requirement that a signal 
be detected in both the Livingston and at least one of the Hanford interferometers 
to construct a triggered search pipeline. 



5.8.1 Selection of Data Quality Cuts for S2 

Playground data from each of the three interferometers was analyzed separately pro- 
ducing a list of inspiral triggers from each interferometer. Only the mandatory data 
quality cuts were used to exclude time from the science mode segments. Each interfer- 
ometer was filtered separately using template banks particular to that interferometer. 
No coincidence was applied between interferometers; data quality cuts were tested in- 
dependently on the three lists of inspiral triggers produced. Table [3] shows the the 
correlation of inspiral triggers with a particular data quality cut for triggers of dif- 
ferent signal to noise. When selecting the data quality cuts we must be aware of 
three constraints. The first is that the data quality cuts are based on data from the 
gravitational wave channel so it is important to ensure that a data quality cut is not 
triggered by a real signal in the data. For this reason we always use caution when 
selecting a cut base on noise in AS_Q. The second constraint is that we do not wish 
to exclude large amounts of data from the analysis. Finally we base our choice on 
advice from the experimental team. A member of the experimental team may decide 
that a cut should be used, even if it does not correlate with inspiral triggers, as any 
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detection made in such a time could not be trusted. Table H] shows the final choice 
of discretionary data quahty cuts and the reasons for them. 

5.8.2 A triggered search pipeline 

During the S2 run, the three LIGO detectors had substantially different sensitivities, 



as can be seen from figure 33 The Livingston interferometer is more sensitive than 
either of the Hanford interferometers. We use this and the requirement that a signal 
be detected in both the Livingston and at least one of the Hanford interferometers to 



construct a triggered search pipeline, summarized in Fig. 34 We search for inspiral 
triggers in the most sensitive interferometer (LI), and only when a trigger is found in 
this interferometer do we search for a coincident trigger in the less sensitive interfer- 
ometers. This approach reduces the computational power necessary to perform the 
search. 

The power spectral density (PSD) of the noise in the Livingston detector is esti- 
mated independently for each LI chunk that is coincident with operation of a Hanford 
detector (denoted LI fl (HI U H2)). The PSD is used to lay out a template bank 
for filtering that chunk, according to the parameters for mass ranges and minimal 
match[72j. The data from the LI interferometer for the chunk is then filtered, using 
that bank, with a signal-to-noise threshold Pl and veto threshold xIl produce 



a list of triggers as described in section 5.3 For each chunk in the Hanford inter- 



ferometers, a triggered bank is created by adding a template if it produced at least 
one trigger in LI during the time of the Hanford chunk. This is used to filter the 
data from the Hanford interferometers with signal-to-noise and thresholds spe- 
cific to the interferometer, giving a total of six thresholds that may be tuned. For 
times when only the H2 interferometer is operating in coincidence with LI (denoted 
LI n (H2 — HI)) the triggered bank is used to filter the H2 chunks that overlap with 
LI data; these triggers are used to test for L1-H2 double coincidence. All HI data 
that overlaps with LI data (denoted LI fl HI) is filtered using the triggered bank for 
that chunk. For HI triggers produced during times when all three interferometers 
are operating, a second triggered bank is produced for each H2 chunk by adding a 
template if it produced at least one trigger found in coincidence in LI and HI during 
the time of the H2 chunk and the H2 chunk is filtered with this bank. These triggers 
are used to search for triple coincident triggers in H2. The remaining triggers from 
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HI when H2 is not available are used to search for Ll-Hl double coincident triggers. 
5.8.3 A directed acyclic graph (DAG) for the S2 pipeline 



In this section we demonstrate how the S2 pipeline in figure 34 can be abstracted into 
a DAG to execute the analysis. We illustrate the construction of the DAG with the 
short list of science segments shown in table [5] For simplicity, we only describe the 
construction of the DAG for zero time lag data. The DAG we construct filters more 
than the absolute minimum amount of data needed to cover all the double and triple 
coincident data, but since we were not computationally limited during S2, we chose 
simplicity over the maximum amount of optimization that could have used. 

A DAG consists of nodes and edges. The nodes are the programs which are 
executed to perform the inspiral search pipeline. In the S2 pipeline, the possible 
nodes of the DAG are: 

1. DATAFIND locates data for a specified time interval on the compute cluster and 
creates a file containing the paths to the input data that other programs can 
read. 

2. TMPLTBANK generates an average power spectral density for a chunk and com- 
putes a template bank for a given region of mass parameter space and minimal 
match. 

3. INSPIRAL filters an analysis chunk using a template bank and generates inspiral 
triggers for further analysis. 

4. TRIGTOTMPLT generated a triggered template bank from the output of the 
inspiral code. 

5. INCA (INspiral Coincidence Analysis) implements the coincidence analysis de- 



scribed in section 5.4 and figure 31 



The edges in the DAG define the relations between programs; these relations are 
determined in terms of parents and children, hence the directed nature of the DAG. 
A node in the DAG will not be executed until all of its parents have been successfully 
executed. There is no limit to the number of parents a node can have; it may be zero 
or many. In order for the DAG to be acyclic, no node can be a child of any node that 
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depends on the execution of that node. By definition, there must be at least one node 
in the DAG with no parents. This node is executed first, followed by any other nodes 
who have no parents or whose parents have previously executed. The construction 
of a DAG allows us to ensure that inspiral triggers for two interferometers have been 
generated before looking for coincidence between the triggers, for example. 

The S2 DAG is generated by a program called the pipeline script, which is an 
implementation of the logic of the S2 pipeline in the Python programming language. 
The pipeline script takes as input the list of science segments for each interferometer, 
with data quality cuts applied. The script reads in all science segments longer than 
2048 seconds and divides them into master analysis chunks. If there is data at the 
end of a science segment that is shorter than 2048 seconds, the chunk is overlapped 
with the previous one, so that the chunk ends at the end of the science segment. An 
option named trig-start-time is set and passed to the inspiral code. No triggers 
are generated before this time and so no triggers are duplicated between chunks. For 
example, the first LI science segment in the fake segment list in table [5] starts at GPS 
time 730000000 and ends at GPS time 730010000. It is divided into the following 



master chunks: 










<AnalysisChunk : 


start 


730000000, 


end 


730002048> 


<AnalysisChunk : 


start 


730001920, 


end 


730003968> 


<AnalysisChunk : 


start 


730003840, 


end 


730005888> 


<AnalysisChunk : 


start 


730005760, 


end 


730007808> 


<AnalysisChunk : 


start 


730007680, 


end 


730009728> 


<AnalysisChunk : 


start 


730007952, 


end 


730010000, trig.start 730009664> 



Although the script generates all the master chunks for a given interferometer, not all 
of them will be filtered. Only those that overlap with double or triple coincident data 
are used for analysis. The master analysis chunks are constructed for LI, HI and H2 
separately by reading in the three science segment files. The full list of master chunks 
for the fake segments is written to a log file. 

The pipeline script next computes the disjoint regions of double and triple coinci- 
dent data to be searched for triggers. 64 seconds is subtracted from the start and end 
of each science segment (since this data is not searched for triggers) and the script 
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performs the correct intersection and unions of the science segments from each in- 
terferometer to generate the following segments containing the times of science mode 
data to search: 

Writing 2 Ll/Hl double coincident segments 
<ScienceSegment : start 730007936, end 730009936, 
<ScienceSegment : start 731001064, end 731002436, 
total time 3372 seconds 



dur 2000> 
dur 1372> 



Writing 2 L1/H2 double coincident segments 
<ScienceSegment : start 730002564, end 730004064, dur 1500> 
<ScienceSegment : start 731004564, end 731005936, dur 1372> 
total time 2872 seconds 



Writing 2 L1/H1/H2 triple coincident segments 
<ScienceSegment : start 730004064, end 730007936, dur 3872> 
<ScienceSegment : start 732000064, end 732002936, dur 2872> 
total time 6744 seconds 

The GPS start and end times are given for each segment to be searched for triggers. 
The script uses this hst of science data to decide which master analysis chunks need 
to be filtered. All LI master chunks that overlap with HI or H2 science data to 
be searched are filtered. An LI template bank is generated for each master chunk 
and the LI data is filtered using this bank. This produces two intermediate data 
products for each master chunk, which are stored as XML data. The intermediate 
data products are the template bank file, Ll-TMPLTBANK-730000000-2048.xml, and 
the inspiral trigger file, Ll-INSPIRAL-730000000-2048.xml. The GPS time in the 
filename corresponds to the start time of the master chunk filtered and the number 
before the . xml file extension is the length of the master chunk. 

All H2 master chunks that overlap with the L1/H2 double coincident data to 
filter are then analyzed. For each H2 master chunk, a triggered template bank is 
generated from LI triggers between the start and end time of the H2 master chunk. 
The triggered bank file generated is called H2-TRIGBANK_Ll-730002500-2048 . xml, 
where the GPS time corresponds to start time of the master H2 chunk to filter. All 
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LI master chunks that overlap with the H2 master chunk are used as input to the 
triggered bank generation to ensure that all necessary templates are filtered. The H2 
master chunks are filtered using the triggered template bank for that master chunk 
to produce H2 triggers in files named H2-INSPIRAL_Ll-730002500-2048.xml. The 
GPS start time in the file name is the start time of the H2 master chunk. 

All HI master chunks that overlap with either the Ll/Hl double coincident data or 
the L1/H1/H2 triple coincident data are filtered. The bank and trigger generation is 
similar to the LI /H2 double coincident case. The triggered template bank is stored in 
a file names Hl-TRIGBANK_Ll-730004000-2048 . xml and the triggers in a file named 
Hl-INSPIRAL_Ll-730004000-2048.xinl where the GPS time in the file name is the 
GPS start time of the HI master chunk. The H2 master chunks that overlap with 
the L1/H1/H2 triple coincident data are described below. 

For each Ll/Hl double coincident segments to search, an inca process is nm to 
perform the coincidence test. The input to inca is all LI and HI master chunks that 
overlap the segment to search. The GPS start and stop times passed to inca are the 
start and stop times of the double coincident segment to search. The output is a file 
names Hl-INCA_LlHl-730007936-2000.xml. The GPS start time in the file name 
is the start time of the double coincident segment. A similar procedure is followed 
for each L1/H2 double coincident segment to search. The output files from inca are 
names H2-INCA_L1H2-731004564-I372.xml, and so on. 

For each L1/H1/H2 triple coincident segment, an inca process is run to create the 
Ll/Hl coincident triggers for this segment. The input files are all LI and HI master 
chunks that overlap with the segment. The start and end times to inca are the start 
and end times of the segment. This creates a file named 

Hl-INCA_LlHlT-730004064-3872 . xml 

where the GPS start time and duration in the file name are those of the triple coinci- 
dent segment to search. For coincidence between LI and an LHO interferometer, we 
only check for time, rft, and mass, dm, coincidence. The parameter k in the effective 
distance cut is set to 1000, so the amplitude cut is disabled. 

For each H2 master chunk that overlaps with triple coincident data, a triggered 
template bank is generated. The input file to the triggered bank generation is 
the inca file for the segment to filter that contains the master chunk. The start 
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and end times of the triggered bank generation are the start and end times of the 
master chunk. This creates a file called H2-TRIGBANK_LlHl-730004420-2048.xinl. 
The H2 master chunk is filtered through the inspiral code to produce a trigger file 
H2-INSPIRAL_LlHl-730004420-2048 . xml. 

For each triple coincident segment to filter, and inca is run between the HI triggers 
from the LIHIT inca and the H2 triggers produced by the inspiral code. The input 
files are the HI inca file Hl-INCA_LlHlT-730004064-3872.xinl and all H2 master 
chunk inspiral files that overlap with this interval. The coincidence is performed as 
follows: 

1. For each HI trigger compute the effective distance of the trigger minus k times 
the effective distance (this is the lower bound on the error allowed in effective 
distance) . 

2. Compute the maximum range of H2 for the trigger mass. 

3. If the lower bound on the HI trigger is further away than can be seen in H2, 
keep the trigger. 

4. If the lower bound on the effective distance of the HI trigger is less than the 
range of H2, but the upper bound is greater, keep the trigger in HI. If a H2 
trigger is found within the interval, store it as well. 

5. If upper bound on the distance of the HI trigger is less than the range of H2, 
check for coincidence. A coincidence check is performed in 5t, 5m, e and k. If 
there is no coincident trigger discard the HI trigger. 

This coincidence step creates two files 

Hl-INCA_LlHlH2-730004064-3872 . xml 
H2-INCA_LlHlH2-730004064-3872 . xml 

where the GPS start time and duration of the files are the start and duration of 
the triple coincident segment. The Ll/Hl coincidence step is then executed again to 
discard any LI triggers coincident with a HI triggered that has been discard by H2. 
The input to the inca are the files 
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Ll-INCA_LlHlT-730004064-3872 . xml 
Hl-INCA_LlHlH2-730004064-3872 . xml 

and the output is the files 

Ll-INCA_LlHlH2-730004064-3872 . xml 
Hl-INCA_LlHlH2-730004064-3872 . xml . 

The HI input file is overwritten as it is identical to the HI output file. Finally, we 
obtain the data products of the search which contain the candidate trigger found by 
the S2 pipeline in these fake segments. The for the fake segments described here, the 
final output files will be: 

Double Coincident Ll/Hl Data 

Ll-INCA_LlHl-730007936-2000 . xml L1-INCA_L1H1-731001064-1372 . xml 
Hl-INCA_LlHl-730007936-2000.xml Hl-INCA_L1H1-731001064-I372.xml 

Double Coincident L1/H2 Data 

Ll-INCA_LlH2-730002564-1500 . xml L1-INCA_L1H2-731004564-1372 . xml 
H2-INCA_LlH2-730002564-1500 . xml H2-INCA_L1H2-731004564-1372 . xml 

Triple Coincident L1/H1/H2 Data 

Ll-INCA_LlHlH2-730004064-3872 . xml Ll-INCA_LlHlH2-732000064-2872 . xml 
Hl-INCA_LlHlH2-730004064-3872 . xml Hl-INCA_LlHlH2-732000064-2872 . xml 
H2-INCA_LlHlH2-730004064-3872 . xml H2-INCA_LlHlH2-732000064-2872 . xml 



As the size of the input science segment files increase, so the number of nodes 
and vertices in the DAG increases, however the algorithm for generating the DAG 
remains the same. 
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Figure 25 : A simple pipeline used to search data from two interferometers for inspiral signals. 
Raw data from interferometer, labeled 1 and 2, is recorded at the observatories. Data quality cuts 
are then apphed to the raw data to discard times when the interferometer was not in a stable 
operating mode. Power spectra generated from the data are used to generate a template bank for 
the inspiral population being searched for. The template banks and interferometer data are used to 
generate inspiral triggers for each interferometer. The triggers for each interferometer are tested for 
coincidence, as a true inspiral signal should be present in both interferometers at the same time, up 
to the time it takes a gravitational wave to travel between the observatories. Other coincidence tests, 
such as waveform parameter consistency, can be applied at this stage. Transient noise sources may be 
detected in auxiliary interferometer channels, for example seismometers. Such channels may be used 
to veto triggers that survive the coincidence analysis but axe coincident with the signature of noise 
in the auxihary channel. Finally we obtain a sample of candidate events for further investigation. 
Each step of the pipeline has many parameters that can be tuned to minimize the false alarm and 
false dismissal rates. 
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Figure 26 : Left Panels: The largest signal-to-noise ratio candidate event seen during the 

search of the LIGO SI data. The top panel shows the signal-to-noise time series, p{t). Notice 
that p{t) is greater than the SI threshold of 6.5 many times in a ~ 5 second interval around 
the candidate event. The center panel shows x^/ {p + S^P^) as a function of time for the val- 
ues of (5^ = 0.03 and p = 8 used in SI. Notice x^/{P + 0.03/9^) is greater than the thresh- 
old of 5 for ~ 5 seconds around the candidate event, but drops below this threshold right at 
the time of maximum p. The inset shows this more clearly for ±0.1 second around the event 
where the threshold is indicated by a dot-dashed horizontal line. The bottom panel shows 
the time series for this candidate event after applying a high-pass filter with a knee frequency 
of 200 Hz. Notice the bursting behavior which does not look like an inspiral chirp signal. 
Right Panels: A simulated injection into the LI data. This example was chosen for comparison 
with the largest signal-to-noise ratio event shown in the left panels since it similar in mass param- 
eters, detected signal to noise and x^- The instrument was behaving well at the time around the 
simulated injection. The top panel shows that p{t) < 6.5 except in close proximity to the signal 
detection time. The center panel shows x^/(p+ O.OBp^) as a function of time. Notice that it is much 
closer to threshold at all times around the simulated injection; this contrasts dramatically with the 
case of the candidate event shown in the left panels. The inset shows this more clearly for ±0.1 
seconds around the injection. The bottom panel shows the time series for this simulated injection 
after applying a high-pass filter with a knee frequency of 200 Hz. The inspiral chirp signal is not 
visible in the noisy detector output. 
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Mandatory Data Quality Cut 



Description 



0UTSIDE_S2 

MISSING_RAW 

DAQ_DROPOUT 

MISSING RDS 

INVALID.TIMING 

CALIB_LINE_NO_RDS_V03 

DAQ_REBOOT 

INVALID CALIB LINE 

NO.CALIB 

LOW.CALIB 



Data is outside of official S2 time interval 

Raw data is missing 

Dropout in data acquisition system 

Data is unavailable for analysis 

Timing believed to be unreliable 

Problem with accessing data for calibration 

One or more data acquisition system rebooted 

Problem with calibration line strength 

Calibration fine turned off 

Calibration fine strength too low 



Discretionary Data Quality Cut 
MICH FILT 

AS_PD_SATURATION 
ASQ_LARGEP2P 

NONSTAND CTRLS 

ASQ_OUTLIER_CLUSTER 

ASQ OUTLIER CORRELATED 

ASQ_LOWBAND .OUTLIER 
ASQ UPPERBAND OUTLIER 



Description 

One or more Michelson control loop 
filters was not in its nominal state 
Antisymmetric port photodiode saturated 
Large peak-to-peak range seen in AS_Q 
at end of lock 

Non-standard controls affecting calibration 
and couplings 

Cluster of large AS_Q outliers 
in short time interval 
Large ASQ outhers correlated with 
outliers in auxiliary IFO channel 
High noise below 100 Hz in AS_Q 
High noise in 100-7000 Hz in AS Q 



Table 2 : Data quality cuts available during the S2 science run and their meanings. Some data 
quality flags monitor human error in the operation of the instrument that make the data unsuitable 

for analysis, such as NO_CALIB and NONSTAND _CTRLS. Others cuts identify hardware or soft- 
ware failures in the operation of the instrument, for example MISSING_RAW and DAQJREBOOT. 
Additional cuts like AS_PD_SATURATION and ASQ_UPPERBAND .OUTLIER monitor the grav- 
itational wave channel, AS_Q, for unusable data. 
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Figure 27 : The loss in signal-to-noise ratio for a 0.2, 0.2 black hole MACHO binary due to 
resampling. The inspiral waveform is generated at 16 384 Hz and injected into data with a typical 
S2 noise curve. The end time of the waveform is at 72.74525 seconds, shown by the vertical line in 
both plots. The top panned shows the signal-to-noise ratio for the data segment when using data 
at the full bandwidth and data resampled to 4096 Hz. The bottom panel shows the same data close 
to the end time of the injection. The loss in signal-to-noise ratio for the inspiral trigger generated 
at 4096 Hz is 0.1%. 
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128 s overlap 



Figure 28 : The algorithm used to divide science segments into data analysis segments. Science 
segments are divided into 2048 s chunks overlapped by 128 s. (Science segments shorter than 2048 s 
are ignored.) An additional chunk with a larger overlap is added to cover any remaining data at 
the end of a science segment. Each chunk is divided into 15 analysis segments of length 256 s for 
filtering. The first and last 64 s of each analysis segment is ignored, so the segments overlap by 
128 s. Areas shaded black are filtered for triggers by the search pipeline. The gray area in the last 
chunk of the science segment is not searched for triggers as this time is covered by the preceding 
chunk, however this data is used in the PSD estimate for the final chunk. 
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Figure 29 : The template bank required to cover the binary black hole MACHO parameter space 
from 0.2 to 1 Mq at a minimal match of 95%. The template bank is generated from the average 
power spectral density of a typical S2 analysis chunk (starting at GPS time 734256712.) The large 
number of templates in the BBHMACHO bank is due to the lager number of cycles of the MACHO 
templates in the sensitive band of the interferometer. The placement of templates outside the mass 
parameter space is required to ensure that any signal that lies in the space has a minimal match 
> 0.95. 
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Figure 30 : The ratio of the known effective distance of an injected signal in the Hanford Observatory 
(LHO) to the known effective distance of an injected signal in the Livingston Observatory (LLO) 
as a function of Greenwich Mean Sidereal Time. The slight misalignment of the interferometers at 
the two different observatories due to the curvature of the earth causes the antenna pattern of the 
detectors to differ. As a result the distance at which a binary system appears is different in each 
detector, even in the absence of noise. The ratio of (^fF(x;tive distances can be significant, so this 
precludes the use of an amplitude cut when testing for inspiral trigger coincidence between different 
observatories. 
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Figure 31 : The test to decide if a trigger in the first detector has a coincident trigger in the second 
detector. If detectors are at different sites, time and mass coincidence are demanded. The effective 
distance cut is disabled by setting k = 1000. If the detectors are at the same site, we ask if the 
maximum distance to which H2 can see at the signal-to-noise threshold pj^2 is greater than the 
distance of the HI trigger, allowing for errors in the measurement of the trigger distance. If this is 
the case, we demand time, mass and effective distance coincidence. If distance to which H2 can see 
overlaps the error in measured distance of the HI trigger, we search for a trigger in H2, but always 
keep the HI trigger even if no coincident trigger is found. If the minimum of the error in measured 
distance of the HI trigger is greater than the maximum distance to which H2 can detect a trigger 
we keep the HI trigger without searching for coincidence. 
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Figure 32 : An auxiliary channel veto investigation of two candidate triggers. Panel (i) shows the 
gravitational wave channel, L1:LSC-AS_Q, high passes above 100 Hz for a an inspiral trigger at GPS 
time 730592784 with a signal-to-noise ratio p = 10.6. The vertical line shows the time of the inspiral 
trigger. Panels (iii) shows the auxiliary interferometer channel L1:LSC-P0B J high passed above 70 
Hz. The vertical lines with circles show the location of glitchMon triggers produced by the noise in 
L1:LSC-P0B J. By excluding inspiral triggers within a time window of these glitchMon triggers, we 
can reduce the number of false event candidates in the pipeline. For contrast, panel (ii) shows an 
gravitational wave channel high passed above 100 Hz for an inspiral trigger at GPS time 734153360. 
This trigger has a similar signal-to-noise ratio, p = 10.9, as the trigger in panel (i). For this trigger 
there does not seem to be a correlated noise source in the auxiliary channel L1:LSC-P0BJ shown 
high passed above 70 Hz in panel (iv). 
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Figure 33 : Typical sensitivities of the three LIGO interferometers during the second LIGO science 
run shown as strain amplitude spectral density, /i/VHz. The smooth solid curve shows the design 
sensitivity (SRD Goal) of the 4 km interferometers and the smooth dashed curve shows the design 
sensitivity of the 2 km interferometer. 
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11400 


11400 


4443 


3922 


3185 



H2 Data Quality Cut 


Ttotal 


T 

-tplay 


T 

done 


p > 8 


p > 10 


p > 12 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


AS_PD_SATURATION 


4 

















MICH_FILT 


64368 


6570 


5648 


1294 


164 


7 



LI Data Quality Cut 


Ttotal 


^play 


-^done 


p > 8 


p > 10 


p>l 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


ASQ_LARGEP2P 


2699 


380 














ASQ_0UTL1ER_C0RRELATED 


840 


60 


60 











AS_PD_SATURAT10N 


646 


61 


10 


813 


119 


6 


M1CH_F1LT 


203539 


21696 


17794 


6393 


497 


32 


NONSTAND.CTRLS 


4020 


843 


18 












Table 3 : The table shows the inspiral triggers generate from science mode data with the mandatory 
data quality cuts applied. For each discretionary data quality cut applied a given interferometer (1), 
the amount of time that would be excluded from the total science mode data by the cut is given 
(2). Since we tune data quality cuts on playground data, the amount of playground time excluded 
is also shown (3) and the amount of playground data analyzed for triggers (4). These may differ for 
reasons explained in section 5.3.2[ The number inspiral triggers generated when a particular data 
quality cut is active is shown for different signal-to-noise thresholds (5-7). To generate the triggers, 
interferometer data was high passed above 50 Hz in the time domain and a low frequency cutoff of 
70 Hz was applied to frequency domain. Template banks were generated with a minimal match of 
0.97 and the signal-to-noise threshold for the matched filter was set to = 8. A veto with 8 
bins applied with a threshold of < 20(8 -I- 0.03^p^). 
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Discretionary Data Quality Cut Applied 



MICH_FILT No 

The cut would exclude a large number of triggers, but would reduce the amount of data in the search 
significantly. It was decided not to apply this cut and to try and exclude false triggers from these 
times by a combination of coincidence, vetoes and reducing the threshold. 

AS_PD_SATURATION Yes 

Clear correlation with inspiral triggers with large signal-to-noise ratios in LI and the study described 



in section 5.2.1 suggest that this should be used. The lack of correlated trigger in HI was due to 



the fact that the playground did not sample any times with photodiode situations. 

ASQ_LARGEP2P No 

A loud inspiral signal could trigger this cut, so it is unsafe for use. 

NONSTAND_CTRLS Yes 

Advice from experimental team advised that detections made during this time could not be trusted. 

ASQ_OUTLIER_CLUSTER No 

Not well correlated with inspiral triggers. 

ASQ_OUTLIER_CORRELATED No 

Not well correlated with inspiral triggers. 

ASQ_LOWBAND_OUTLIER No 

Not well correlated with inspiral triggers. 



ASQ_UPPERBAND_OUTLIER Yes 

Times with high upper band noise in HI are clearly correlated with high signal-to-noise ratio triggers. 
In order to prevent the cut from begin triggered by real signals we also require that the cut is on 
for more that 180 seconds. The longest inspiral signal in the S2 analysis is 52 seconds. 

Table 4 : The final selection and justification of discretionary data quality cuts for the S2 binary 
neutron star and binary black hole MACHO searches. 



140 



L1 GW Data 
L1 n (H1 u H2) 



H1 GW Data 




Template Bank 
Generation 



Matched Filter 
p2, x2 threshold 



Triggers to 
Template Bank 



Matched Filter 
p2, x2 threshold 



H2 GW Data 



Triggers to 
Template Bank 



L1/H1 Trigger 
Coincidence 



H2 GW Data 
L1 n H1 n H2 



Triggers to 
Template Bank 



Matched Filter 
p2, x2 threshold 



H1/H2 Trigger 
Coincidence 



+ 



L1/H1 
Triggers 



L1/H1/H2 
Triggers 



Matched Filter 
p2, x2 threshold 



L1/H2 Trigger 
Coincidence 



L1 Bean 
Pickoff 


1 Splitter 
Channel 






Glitch Monitor 






Veto Times 



L1/H2 
Triggers 



Figure 34 : The inspiral analysis pipeline used to determine the reported upper limit. Lin (H1UH2) 
indicates times when the LI interferometer was operating in coincidence with one or both of the 
Hanford interferometers. LI D HI indicates times when the LI interferometer was operating in 
coincidence with the HI interferometer. LI D (H2 — HI) indicates times when the LI interferometer 
was operating in coincidence with only the H2 interferometer. The outputs of the search pipeline 
are triggers that belong to one of the two double coincident data sets or to the triple coincident data 
set. 
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Interferometer Start End Duration 

LI 730000000 730010000 10000 

LI 731001000 731006000 5000 

LI 732000000 732003000 3000 

HI 730004000 730013000 8000 

HI 731000000 731002500 2500 

HI 732000000 732003000 3000 

H2 730002500 730008000 5500 

H2 731004500 731007500 2500 

H2 732000000 732003000 3000 



Table 5 : The fake science segments used to construct the DAG shown in figure 
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Chapter 6 



Hardware Signal Injections 



Gravitational radiation incident on the LIGO interferometers from an inspiralling 
binary will cause the test masses to move relative to each other. This produces a 



differential change in length of the arms as described in section 2A_ Injection is 
the process of adding a waveform to interferometer data to simulate the presence of 
a signal in the noise. We use injections to measure the performance of the binary 



inspiral analysis pipeline as described in section |5.7[ Software injections, which add 
a simulated signal to the data after it has been recorded, are used for efficiency 
measurements. Since they performed a posteriori the interferometer is not affected 
while it is recording data. Alternatively, a simulated signal can be added to the 
interferometer control system to make the instrument behave as if an inspiral signal 
is present. The interferometer Length Sensing and Control system has excitation 
points which allow arbitrary signals to be added into the servo control loops or to the 
drives that control the motion of the mirrors[74]. We call this hardware injection; the 



data recorded from the instrument contains the simulated signal. Figure [36] shows the 
hardware injection points on a schematic diagram of the interferometer and length 
sensing and control loop. 

Analysis of hardware injections allows us to ensure that the analysis pipeline is 
sensitive to real inspiral signals and validates the software injections used to test the 
pipeline efficiency. In order to perform an accurate upper limit analysis for binary 
inspirals, we must measure the efficiency of our pipeline. That is, we inject a known 
number of signals into the pipeline and determine the fraction of these detected. 
Injecting signals into the interferometer for the duration of a run is not practical 
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and would contaminate the data, so we use the analysis software to inject inspiral 
signals into the data. By comparing software and hardware injections we confirm that 
software injections are adequate to measure the efficiency of the upper limit pipeline. 

Hardware injections provide a very complete method of testing the inspiral detec- 
tion pipeline. By recovering the physical parameters of an injected signal, we test our 
understanding of all aspects of the pipeline, including the instrumental calibration, 
the filtering algorithm and veto safety. We injected inspiral signals immediately after 
the first LIGO science run (SI) in September 2002. The resulting data was analyzed 
using the software tools used to search for real signals. In this chapter, we describe 
the results of analysis of the SI hardware injections. The analysis pipeline used in SI 
differs from that used in S2[T2]. Here we are examining the response of the filtering 
code to the hardware injections, however, and so the differences between the SI and 
S2 pipelines are unimportant. 



6.1 Injection of the Inspiral Signals 

To inject the signals, we generate the interferometer strain h{t) produced by an in- 
spiralling binary using the restricted second order post-Newtonian approximation in 
the time domain[2lj. The LSC calibration group supplies a transfer function T(/) 
which allows us to construct a signal g{t) that produces the desired strain when it is 
injected into the interferometer. The transfer function T(/) should be identical to the 



actuation function A(f) described in section 2.2.3, however in SI this was simplified 



to contain only the pendulum response of the mirrors, given by 



LP 



nf) = (6-1) 



where L is the length of the interferometer, C is the calibration of the excitation 
point in nm/count and /o is the pendulum frequency of the test mass. Damping 
is neglected as it is unimportant in the LIGO frequency band. The code used to 
generate the hardware injections is the same as that used for software injections; only 
the transfer function used to generate the injected signal differs since we are injecting 
into the control signal g rather than the error signal v. 

During SI, we injected signals corresponding to an optimally oriented binary. 
Injections of a 1.4 Mq inspiralling binary at distances from 10 kpc to 80 kpc were 
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used to test the neutron star analysis. We also injected signals from a 1.4, 4.0 Mq 
binary and several 1.4, 1.4 binaries at closer distances. These signals were injected 
into the differential mode servo and directly into an end test mass drive. 



6.2 Detection of the Injected Signals 



Figure [37] shows the events generated by processing 4000 seconds of data from the Liv- 
ingston 4 km interferometer (LI) on 10 September 2002 during the post-run hardware 
injections. 

The first set of injections were large amplitude signals used to verify the inspirals 
were being correctly injected. We ignore these and concentrate on the second set, 
which were at more appropriate distances. We only consider the 1.4 solar mass 
inspiral injections, as the 1.4,4.0 injection lies outside the template bank space used 
in the SI binary neutron star analysis. 

It can be seen that all of the hardware injections are identified as candidate events 
since they have high signal-to- noise ratios and values of the test lower than 5, which 
was the threshold used in the SI analysis pipeline [T^. Some of the 1.4,4.0 injections 
are also flagged for further investigation as they cause templates inside the bank to 
ring, but have high values as they are not exactly matched. 

Since we know the exact coalescence time of the injected waveform, we can com- 
pare this with the value reported by the search code and ensure that the search code 
is reporting the correct time. The known and measured parameters for the second 
set of 1.4, 1.4 Mq injections are shown in table |6} The raw data is resampled to 4096 
Hz before being filtered. For each of the signals injected, we were able to detect the 
coalescence time of the injection to within one sample point of the correct value at 
4096 Hz, which is consistent with the expected statistical error and confirms that the 
pipeline has not introduced any distortion of the signals. 



6.3 Checking the Instrumental Calibration 

Calibration measurements of the interferometers were performed before and after the 
run; these are the reference calibrations. In general, the calibration changes due to 
changes in the alignment on time scales of minutes. This variation is encoded in the 
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parameter a which is monitored using a sinusoidal signal injected into the detector 



(see section 2.2.3 ). a is used as input to the data analysis pipeline and varied between 
0.4 and 1.4 during SI. Data in SI was analyzed in 256 second segments. For each 256 
seconds of data starting at time to? we construct the calibration, -R(/; to) by using 
a(to) and a reference calibration. to) is then used to calibrate the 256 seconds 
of data. 



Figure 38 shows a set of injections into the Livingston interferometer analyzed 
with different calibrations generated by varying the value of a. We expect that 
the signal-to-noise varies quadratically and the effective distance varies linearly with 
changes in a [75]. This is confirmed by the injections. There is no single value of a 
that gives the correct effective distance for all the injections; this is consistent with 
the estimated systematic errors in the calibration. Unfortunately the calibration line 
was not present during the time the hardware injections were performed, so we cannot 
directly compare a measured calibration with the result of the injections. 

6.4 Safety of Vetoes 

During construction of the the inspiral pipeline we considered using inspiral triggers 
found in auxiliary interferometer channels as vetoes on triggers in the gravitational 
wave channel. Concern was raised that a real inspiral signal may couple between these 
channels and a real signal may be inadvertently vetoed. To check this, we examined 



coupling between the channels at the time of an injection. Figure [39] shows the power 
spectra of the gravitational wave channel, LSC-AS_Q, and the auxiliary channels that 
we considered using as vetoes during SI: LSC-AS_I, LSC-REFL_I and LSC-REFL_Q. 
The injected inspiral signal can clearly be seen coupling to the auxiliary channel LSC- 
AS_I, but there is no obvious coupling between the injected signal and LSC-REFL_I 
or LSC-REFL_Q. This led us to discard LSC-AS J as a possible veto channel in the SI 
analysis. Similar studies have been performed for the S2 data when auxiliary channels 
are proposed as veto channels. 
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Figure 36 : A schematic diagram of the LIGO interferometer showing the injection points used in 
SI hardware injections. Inspiral signals were injected either directly into the end test mass drive of 
one arm or into the differential mode servo, and this into both arms. Care was taken to ensure that 
the correct transfer function, T{f), was used in each case. 



End time of Injection 


End Time of Detection 


P 




04 : 35 : 12.424928 


04 : 35 : 12.424927 


11.623546 


1.653222 


04 : 36 : 42.424928 


04 : 36 : 42.425171 


20.230101 


1.671016 


04 : 38 : 12.424928 


04 : 38 : 12.424927 


37.488770 


0.443966 


04 : 39 : 42.424928 


04 : 39 : 42.424927 


69.815262 


1.375486 



Table 6 : Hardware injection events found by the inspiral analysis pipeline. End time of injection 
is the known end time of the injected signal and end time of detection is the end time of the signal 
as reported by the analysis pipeline. Times are Universal Time (UTC) on 10 September 2002. The 
values of signal-to- noise ratio p and veto are given for each event. 
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Figure 37 : The candidate events generated by processing 4000 seconds of data from the Livingston 
4 km interferometer through the SI analysis pipehne. This data included two sets of injections; 
the known coalescence times are indicated by the dashed vertical lines. The signal-to-noise ratio is 
plotted and the value of the veto is shown next to the candidate event. 
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value of calibration parameter a 



Figure 38 : Each curve corresponds to a hardware injection at the given GPS time. We re-analyze 
each injection with different caHbrations to show how the detected quantities vary with a. The 
upper plot shows the ratio of signal-to-noise ratio, p, to its maximum value, Pmax- The lower plot 
shows the ratio of the detected distance to the known distance of the hardware injection. 
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Figure 39 : Power spectra of the gravitational wave channel LSC-AS_Q and the auxiliary channels 
LSC-ASJ, LSC-REFLJ and LSC-REFL_Q during a hardware injection. The broad peak in the 
spectrum is the inspiral signal and the two power spectra taken at subsequent times show it sweeping 
across the band as the frequency of the inspiral signal increases with time. 
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Chapter 7 

The Rate of Binary Black Hole 
MACHO Inspirals in the Halo 



In this chapter, we present the resuhs of a search for gravitational waves from the 
inspiral of binary black hole MACHOs in data from the second LIGO science run 
(called S2). The goal of the search is the detection of the gravitational waves. In the 
absence of a detection, however, we place an upper limit on the rate of inspiralling 
BBHMACHOs. This limit may be compared to the predicted rate of 5 x 10^^ x 2^^ 
discussed in chapter |3| 

Analysis of the full S2 data set for gravitational waves from inspiralling binary 
black hole MACHOs is complete and the result of this search will appear in [76] . 
Since this result is currently embargoed pending LIGO Scientific Collaboration in- 
ternal review, we instead present the result of the search on the playground data. 
No gravitational waves from BBHMACHO inspirals were found in the playground, 
so in section |7.5 we compute an upper limit on the rate of binary black hole MA- 



CHO inspirals in the playground data. Although this result is statistically biased, as 
it is computed from data used to tune the pipeline, it allows us to make a reason- 
able prediction of the upper limit available using the full S2 data and assuming no 
BBHMACHO signals are detected in the full data set. 



In section 7.1 we describe the data sample used in the analysis. Section 7.2 



describes how the parameters of the search listed in the previous chapter were tuned 



on the playground data. Section 7.3 described the Monte Carlo simulations used to 



measure the efficiency of the pipeline. In section 7.4 we describe the background 
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observed in the S2 data. 



7.1 The Second LIGO Science Run 

All three LIGO detectors operated during the second science run, referred to as S2, 
which lasted for 59 days (1415 hours) from February 14 to April 14, 2003. Although 
the detectors were manned by operators and scientific monitors around the clock, the 
amount of data flagged for scientific analysis was limited by environmental factors 
(especially high ground motion at LLO and strong winds at LHO), occasional equip- 
ment failures, and periodic special investigations. The total amount of science data 
obtained was 536 hours for LI, 1044 hours for HI, and 822 hours for H2. 

The analysis described in this thesis uses data collected while the LLO detector 
was operating at the same time as one or both of the LHO detectors in order to 
make use of the triggered search pipeline. Science mode data during which both HI 
and H2 were operating but LI was not, amounting to 383 hours, was not used in 
this analysis because of concerns about possible environmentally-induced correlations 
between the data streams of these two co-located detectors. This data set, as well 
as data collected while only one of the LIGO detectors was in science mode, will be 



combined with data from the third LIGO science run in a future analysis. Figure 40 
shows a breakdown by interferometer of the data recorded during S2. The data used 
in this search is indicated by the shaded region. 



7.2 Tuning the Analysis PipeUne 

The entire analysis pipeline was explored first using the playground data set in or- 
der to tune the the various thresholds and other parameters. The goal of tuning 
the pipeline is to maximize the efficiency of the pipeline to detection of gravitational 
waves from binary inspirals without producing an excessive rate of spurious candidate 
events. In the absence of a detection, a pipeline with a high efficiency and low false 
alarm rate allows us to set the best upper limit. It should be noted, however, that 
our primary motivation is to enable reliable detection of gravitational waves. The 
efficiency is measured by Monte Carlo simulations in which signals from the hypo- 
thetical population are added to the data and then sought. This approach accounts 
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for any systematic error associated with the methods used in our pipehne. Note that 
another factor in the tuning of the pipehne are the available computational resources. 
We would like to be able to complete the search in less time than the length of the 
data being analyzed, so that real-time searches are possible when the interferometers 
are taking continuous data. For this reason, certain tuning decisions are based on 
the computational efficiency of the pipeline; these decisions will be clearly identified 
below. 

Prior to commencing the binary black hole MACHO search, a search for inspi- 
ralling binary neutron stars (BNS) was conducted on the S2 data using the pipeline 
described in chapter [5]_14j. The mechanics of the BNS search are very similar to 
those described here, except that the template bank covers binaries with 1.0 Mq < 
mi, m2 < 3.0 Mq, where mi and 1712 are the masses of each object in the binary. Since 
the BNS and binary black hole MACHO searches are very similar, and share the same 
playground data, we may use the parameters of the BNS search (which were tuned 
on the playground data) as a starting point for tuning the binary black hole MACHO 
search. 

There are two sets of parameters that we are able to tune in the pipeline: (i) 
the single interferometer parameters which are used in the matched filter and 
veto to generate inspiral triggers in each interferometer, and (ii) the coincidence 
parameters used to determine if triggers from two interferometers are coincident. 
The single interferometer parameters include the signal-to- noise threshold p*, the 
number of frequency sub-bands in the statistic p, the cut threshold S*, and the 



coefficient on the signal-to-noise dependence of the x cut, i.e. 5 in equation (4.151 ). 
These are tuned on a per-interferometer basis, although some of the values chosen 
are common to two or even three detectors. The coincidence parameters are the time 
coincidence window 6t for triggers, the mass parameter coincidence window 6m and 



the effective distance cut parameters e and k in equation. (5.8). Due to the nature 
of the triggered search pipeline, parameter tuning was carried out in two stages. We 
first tuned the single interferometer parameters for the primary detector (LI). We 
then used the triggered template banks (generated from the LI triggers) to explore 
the single interferometer parameters for the less sensitive Hanford detectors. Finally 
the parameters of the coincidence test were tuned. 
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7.2.1 Template Bank Generation 

Recall that the number of templates needed to cover a given region of parameter space 
at a specified minimal match fluctuates as the shape of the noise power spectrum 
changes. The greater the sensitivity at low frequencies, relative to higher frequencies, 



the larger the template bank (see section 5.3.1). The computational resources avail- 
able for the MACHO search are limited and the computational cost is proportional 
number of templates in the bank. We therefore tuned the template bank parameters 
to allow the search to be completed within the available resources. 

Due to the algorithm used to construct the template bank[72], the smallest mass 
template in the bank will be the equal mass binary (mmin, "^min), where is the 



(user specified) minimum binary component mass. Figure 41 shows the size of the 
template bank necessary to cover each playground analysis chunk at a minimal match 
of 0.97% for several values of mmin- The maximum binary component mass in each 
case is IMq, so the largest mass binary in the bank parameter space is (1,1)M0. 
For fixed mmin the number of templates remains reasonably constant over the course 
of the S2 run, but there is a large variation in the number of templates required as a 
function of lower mass. The scaling of template number as a function of lower mass 
is consistent with that described in |72] . 

As described in chapter [3| the MACHO mass range measured by microlensing is 
0.15 M0 to 0.9 M0 at 95% confidence. It would therefore be desirable for the binary 
black hole MACHO inspiral search to cover a region of mass parameter space slightly 
larger than this, say 0.1 Mq to 1.0 M0. It can be seen, however, that almost an 
order of magnitude more templates are needed to decrease the lower boundary of the 
mass parameter space from 0.2 Mq to 0.1 Mq. Therefore, given the computational 
resources available for the BBHMACHO search, the lowest mass template in the bank 
was set to 0.2 Mq. Note that for the final search, the match of the template bank 
was also lowered to 0.95% to further decrease the number of templates to an average 
of 14179 per analysis chunk over the S2 run. The latter choice is be justified by the 
sensitivity to galactic binary black hole MACHOs in S2, as will be seen below. The 



size of the inspiral template bank used is shown in figure 42 
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7.2.2 Interferometer Sensitivity and Signal-to-noise Thresholds 

The noise power spectrum also determines the sensitivity of the interferometer to 
binary inspirals. We can quantify the sensitivity in terms of the distance to which we 
can see an optimally oriented binary inspiral at a given signal-to-noise ratio. This is 
the maximum distance at which the interferometer can detect a binary (at this signal- 
to-noise ratio), since the gravitational wave strain in the interferometer is a maximum 
when the binary is optimally oriented. The maximum inspiral ranges for an optimally 



oriented binary at signal-to-noise ratio p* = 8 are shown in figure |42j The distance 
to which we can detect an optimally oriented binary is also a function of the mass of 
the binary, scaling as jj^^'^M^^^; figure 42 shows the ranges for a (0.5,0.5) Mq and a 
(0.1, 0.1) Mq binary. 

Notice that there are no times when either of the LHO interferometers are more 
sensitive than the LLO interferometer and so demanding triggers are always present 
in the most sensitive interferometer means that they are required to be found in 
LI. In fact the LLO interferometer has a significantly larger range than either of 
the LHO interferometers, at times being sensitive to BBHMACHO inspirals in An- 
dromeda at around 0.7 Mpc. Since we require coincidence between LI and one of the 
LHO interferometers to make a detection, however, we are restricted to a search for 
BBHMACHOs in the Galactic halo. 



Based on the sensitivity plots shown in figure 42, we set the signal-to-noise thresh- 
old to 7 in all three interferometers; we justify this as follows. All three interferometers 
are sensitive to optimally oriented inspirals with p* > 8 at distances greater than the 
size of the Galactic halo. A binary black hole MACHO in the Galaxy may have an 
unfavorable orientation, however, causing it to appear at a large effective distance. 
For this reason, we want to set the signal-to-noise ratio threshold as low as possible 
without producing an excessive false alarm rate. Lowering the signal-to-noise thresh- 
old has a computational impact on our search: when the signal-to-noise ratio for a 
template crosses threshold, we perform the veto which requires p additional com- 
plex inverse FFTs, where p is the number of frequency bins used in the veto. If we 
set the signal-to-noise threshold too low, we may exceed the available computational 
resources due to the extra operations required to perform the veto. In fact this is 
what happens with the S2 binary black hole MACHO search since the template banks 
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are so large. Whereas in the S2 BNS search we were able to lower the signal-to-noise 
threshold to 6, the BBHMACHO search is limited to a signal-to- noise threshold to 7. 

7.2.3 Tuning the Veto Parameters 

Recall from section 4.8.2 that the veto thresholds on 

x'<xliP + p'5'), (7.1) 

where p is the signal-to-noise ratio of the signal and 5^ is a parameter chosen to be 
reflect the largest expected mismatch that a true signal will have with the templates 
in the bank. The initial parameters used for the veto, based on tuning of the BNS 
search, were p = 15 and 5^ = 0.04 with the threshold set to x^(Ll) = 5.0 in the LI 
interferometer and x^(Hl) = x^(H2) = 12.5 in the LHO interferometers. 



Figure 43 illustrates tuning of the veto on the HI playground triggers. If the 
interferometer noise is Gaussian, then the square of the signal-to-noise ratio should 
be distributed with two degrees of freedom. This means that the histogram of 
the signal-to-noise ratio of the triggers should be a monotonically decreasing function 



of p. We can see from figure |43} however, that there is an excess in the number of 
triggers with p ^ 8.5; this suggests some non-Gaussian behavior in the data that we 
would like to remove. It can be seen that decreasing the threshold xl to 5 removes 
this hump in the distribution and decreases the signal-to-noise ratio of the loudest 
event from p = 13.8 to p = 10.7. This suggests that lowering the threshold is 



desirable. Figure 44, which shows the results of a small Monte Carlo simulation, 
demonstrates the danger of making such decisions without reference to the detection 
efficiency, however. For each simulated signal added to the data, the figure shows 
whether or not it was detected in the HI data using the initial choice of veto 
parameters (5^ = 0.04, x^(Hl) = 12.5). Several injections are missed at effective 
distances well within the range of the HI interferometer. Follow up investigations of 
these missed triggers show that, although they have very large values of signal-to- noise 
p ~ 10^ — 10^, they are missed because they fail the veto. This is caused by the 
mismatch between the signal and the template. The results of this study imply that 



we should loosen the x veto, in contradiction to the results suggested by figure |43 
Notice, however, that the excess of HI triggers occurs at low values of p and the missed 
injections are at higher values of p. Figure 45 shows the values of x^/(p + p"^^"^) ^-^id 
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p for the detected LI injections and the detected HI injections before coincidence 
is apphed. These triggers are taken from the same Monte Carlo simulation as the 



triggers shown in figure 44 The values of j {p + p'^S^) observed for the detected HI 
injections are significantly higher than those for the LI injections. This suggests that 
we should increase the parameter 6^, which has the effect of decreasing the value of 
X^/iP + p"^^"^) foi' triggers with high signal-to-noise ratios. If we increase 5^ we may 
be able to decrease the value of x^/{p + p"^^^) to remove the excess of triggers in HI, 
without adversely affecting the pipeline detection efficiency. After several iterations, 
the values 5^ = 0.4, x*(Ll) = 3.1, X*(H1) = 5.0 and x*(H2) = 10.0 were chosen. 



Figure [46] shows the values of S and p for the detected LI injections and the detected 
HI injections with these new parameters. Notice that the values of x^/{p + P^^^) 
for the loud signal-to-noise triggers are considerably lower than before. It appears 



from the results in figure 46 that it would be possible to reduce x*(Hl) further to 
2.51, without loss of efficiency. No coincident triggers survived the pipeline in the 
playground data with a threshold of x*(Hl) = 5.0, however, so it was decided not to 
reduce this threshold further. 

7.2.4 Coincidence Parameter Tuning 

After the single interferometer parameters had been selected, the coincidence param- 
eters were tuned using the triggers from the single interferometers. As described in 
section [6| the coalescence time of an inspiral signal can be measured to within < 1 ms. 
The light travel time between observatories is 10 ms, so 6t was chosen to be 1 ms for 
LHO-LHO coincidence and 11 ms for LHO-LLO coincidence. The mass coincidence 
parameter was initially chosen to be 6m = 0.03, however testing with the binary 
neutron star search showed that this could be set to 6m = 0.0 {i.e. requiring the 
triggers in each interferometer to be found with the exact same template) without 
loss of efficiency. 

Having tuned the time and mass parameters, we tune the effective distance pa- 
rameters K and e. Initial estimates of e = 2 and k = 0.2 were used for testing, 
however it was discovered that many injections were missed using these thresholds to 
test for LLO-LHO coincidence. This is due to the fact that the detectors are slightly 
misaligned, so the ratio of effective distance of a trigger between the two observa- 



tories can be large for a significant fraction of the population, as shown in Fig. 30 
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As a result, we disabled the effective distance cut for triggers generated at different 
observatories. A study of simulated signals injected into HI and H2 interferometers, 
both located at the LIGO Hanford Observatory, suggested using values of enn = 2 
and = 0.5. Note that, as described above, we demand that an Ll/Hl trigger 
pass the H1/H2 coincidence test if the effective distance of the trigger in HI is within 
the maximum range of the H2 detector at threshold. 



7.3 Results of Injection Monte Carlo Simulations 

The final parameter values chosen are shown in table [7} Once fixed, an injection 
Monte Carlo was performed to measure the efficiency of the search pipeline. For each 
playground interval in the S2 data, an inspiral waveform was generated using the 



population model described in section 3.4, These signals has masses between 0.1 Mg 







and 1.0 Mq. Each inspiral signal was injected into the data at a random time during 
a unique playground interval and the data analyzed through the full pipeline with 
the final set of parameters. Four separate simulation runs were performed giving a 



total of 849 injections in the analyzed playground data. Figure 47 shows the results 
of this simulation in the (mi, 7/12) plane. The figure shows which of the simulated 
signals were detected and which were missed by the full pipeline in the double and 
triple coincident playground data. Also shown in the figure is the effective distance 
at which these signals were injected in the LHO interferometers, since it is generally 
the less sensitive detector that limits the detection efficiency. Although the template 
bank only coves the region above and to the right of the red lines at 0.2 Mq, the 
upper plot shows that some signals are detected with component masses between 
~ 0.15-0.20 Mq. This is a direct result of increasing 6"^ in the veto, equation 



(7.1), allowing loud, but slightly mismatched signals to be detected. The lower plot 



in figure 47 shows the injections that are missed by the pipeline. Injections missed 
from triple, Ll-Hl double and L1-H2 double coincident data are shown with a star, 
an upward triangle and a downward triangle respectively. These missed signals are 
color coded with the injected effective distance in the LHO detectors. Injections are 
only missed when their effective distance is comparable or greater than the ranges 



for the LHO detectors shown in figure 42 there are no anomalous missed injections. 



Figure 48 shows the efficiency of the search as a function of chirp mass. As expected. 
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given the strength of the BBHMACHO signals and the sensitivity of the detectors, 
the efficiency is e ~ 1 for > 0.35 with the small loss of efficiency coming from 
systems with an unfavorable orientation. The efficiency drops for Ai below 0.35 due 
to the combined effect of the signals becoming weaker as the chirp mass decreases and 
falling outside the region of good template bank coverage. Notice that there appears 
to be an anomalous value of e at 7V1 ~ 0.2; this appears to be associated with a 
dearth of injections in this mass range. Large scale Monte Carlos, with many more 
injections, are currently being performed to explore this region of parameters space 
and investigate this anomaly in the full S2 data set. 

A comparison of the inspiral parameters recovered by the pipeline and the known 
injection parameters from signals injected in the Monte Carlo simulation are shown 



for LI, HI and H2 in figures 49, 50 and 51 respectively. It can be seen that the effec- 
tive distance recorded by the search code is unbiased in all cases, and can typically be 
recovered to an accuracy of ~ 10%. This is comparable to the 10% distance uncer- 
tainty to nearby galaxies. Using the measured coalescence phase, effective distance 
and difference in time of arrive at the two detectors, it is possible to gain informa- 
tion about the location of the signals, however a comprehensive study of this has not 
been performed for the S2 data. For all interferometers, the chirp mass is recovered 
extremely well with an accuracy of 0.1%. This consistent with the results quoted in 
[27] and is encouraging for the parameter measurement in the case of a detection. 
Although there appears to be a bias in the measure end time of the signal, this is 
due to the fact that the current implementation of the filtering code measures the 
end time of the template, not the coalescence time of the binary. The injected signals 
are generated with a time domain waveform generator [77] and recovered with the 
stationary phase waveforms described in chapter |4j There is a slight difference in the 
frequency at which these waveforms terminate: the time domain waveforms are ter- 



minated when the post-Newtonian phase evolution of equation (2.110) can no longer 
be evolved and the stationary phase waveforms are terminated when they reach the 
gravitational wave frequency of a test particle in the innermost stable circular orbit 
of Schwarzschild spacetime. As a result of this there is a small mass dependent offset 
between the end times of the waveform as recorded by the injection code and the 
filtering code. This will not affect the time coincidence test, however, as we demand 
that the templates have the same mass parameters in all the detectors. As a result 



160 



this time offset will be identical between detectors and we can apply coincidence. 
Changes to the filtering code are planed to remove this offset in time measurement. 

7.4 Background Estimation 

We estimate the background rate for this search by introducing an artificial time 
offset, or lag, At in the triggers coming from the Livingston detector relative to the 
Hanford detector. The time-lag triggers are then fed into subsequent steps of the 
pipeline. The triggers that emerge from the end of the pipeline are considered a 
single trial representative of an output from a search if no signals are present in the 
data. By choosing a lag of more than 20 ms, we ensure that a true gravitational wave 
will not be coincident in the time-shifted data streams. We do not time-shift the two 
Hanford detectors relative to one another since there may be real correlations due to 
environmental disturbances. If the times of background triggers are not correlated at 
the sites, then the background rate can be measured. A total of 20 time-lags were 
analyzed to estimate the background. Note that the time lags use all the data and 
are not restricted to playground. The resulting distribution of time-lag triggers in 



the (ph,Pl) plane is shown in figure 52; the distribution of background triggers and 



injected signals are compared in figure 53 It can be seen that the signal-to-noise 



ratios of background triggers are higher in the LHO interferometers than in the LI 
interferometer, whereas the signal-to-noise ratios of injections are louder in the LI 
interferometer. This distribution suggested the form of a "coherent" signal-to-noise 
ratio p which gives a factor of 2 more significance to the signal-to-noise ratio in LI 
compared to the signal-to-noise ratio in LI. Based on the studies of the background 
and injections, we chose a combined signal-to-noise statistic 

i 

4 



= pI + (7.2) 



7.5 Upper Limit on BBHMACHO Inspiral in the S2 Play- 
ground Data 



After the data quality cuts, discarding science segments with durations shorter than 
2048 sec, and application of the instrumental veto in LI, a total of 35.2 hours of 
playground data remained; 22.0 hours of triple-detector data, 10.2 hours of Ll-Hl 
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and 3.97 hours of L1-H2. The analysis of the playground data produced no double or 
triple coincident inspiral candidates. 

To determine an upper limit on the event rate we use the loudest event statisticf78\ 
which uses the detection efficiency at the signal-to-noise ratio of the loudest trigger 
surviving the pipeline to determine an upper limit on the rate. If no triggers survive 
the pipeline, we use the signal-to-noise threshold as the loudest event. Suppose the 
population of sources produces Poisson-distributed events with a rate 71 per year per 
Milky Way Equivalent Galaxy (MWEG) and Ng{p*) is the number of MWEGs to 
which the search is sensitive at p > p*. Then the probability of observing an inspiral 
signal with p > p*, given some rate TZ and some observation time T, is 

P{p > p*-n) = l- e-^™G(P*)_ (73) 

A trigger can arise from either an inspiral signal in the data or from background. If 
Pft denotes the probability that all background triggers have signal-to-noise ratio less 
than p*, then the probability of observing either an inspiral signal or a background 
trigger with p > p* is given by 

P{p > p*;TZ,b) = l- Pi,e-'^™'^^P*\ (7.4) 

Given the probability Pb, the total observation time T, and the number of Milky Way 
equivalent galaxies Nq to which the search is sensitive, we find that the rate of binary 
black hole MACHO inspirals per MWEG is 

2.303 + In n 

with 90% confidence. This is a frequentist upper limit on the rate. For 71 > 7Zqo%, 
there is more than 90% probability that at least one event would be observed with 
SNR greater than pmax- 

Since no coincident events were observed in the playground data, we determine 
the rate by measuring the efficiency of the pipeline at the (combined) signal-to-noise 
threshold ^ 

pLx = 7' + J = 61-25. (7.6) 

This is a conservative limit on the rate. If we lowered the signal-to-noise thresholds 
until we observed the loudest trigger, then this trigger will have a signal-to-noise ratio 
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less than the value given in equation (7.6); the measured efficiency will be greater 
than or equal to that used here so the rate will be less than or equal to which we 
quoted here. Furthermore, since we only have a small number of time lags, we neglect 
the background term Ph. Dropping this term will give a conservative value for upper 
limit [78] . 

If we restrict the mass parameters of the injected signals to those with a component 
in the range 0.15 to 1.0 M©, the mass range suggested by microlensing observations, 
we find that efficiency of the search pipeline at p* = -\/61.25 is 

e (p* = ^61:25) = ^ = 0.915. (7.7) 

The observation time for the playground data is T = 35.2 hours = 4 x 10^'^ yr and so 
the upper limit on the rate of binary black hole MACHO inspirals in the playground 
data is 

2 303 

= 0,915 X 4 X 10-^ = f^-**' 
The amount of non-playground data in the full S2 data, again discarding science 

segments with durations shorter than 2048 sec, and application of the instrumental 
veto in LI, there is a total of 345 hours of non-playground data; 225 hours of triple- 
detector data, 90 hours of Ll-Hl and 30 hours of L1-H2. Assuming that the signal- 
to-noise ratio of the loudest event in the full data is comparable to = 61.25, which 
is suggested by the background triggers, we may estimate the achievable upper limit 
as 

2 303 

= 0,915X3.9X10-^ = P-^' 
This upper limit is three orders of magnitude larger than the upper bound on the 

rate of i? = 0.05 x 2^^ yr~^ MWEG~^, so we are not yet able to constrain the fraction 

of the Galactic halo in BBHMACHOs. The sensitivity of the interferometers during 

the S2 search is roughly an order of magnitude from design sensitivity, however. For 

a (0.5, 0.5) Mq binary, at design sensitivity the interferometers will be sensitive to 

~ 50 MWEG, so in one year of data taking, assuming no detections have been made, 

an upper limit on rate of BBHMACHO inspirals oi R = 3 x 10"^ yr"^ MWEG"^ 

are possible (assuming Ph = 0.5) which will significantly impact the theoretical rate 

estimates. Additionally, since the range of the search scales as a function of the mass, 

the possibility for detecting more massive BBHMACHO binaries increases as more 

galaxies become accessible and the rate can be further constrained. 
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Figure 40 : The Venn diagram shows the number of hours that each detector combination was 
operational during the S2 run. The upper number gives the amount of time the specific instruments 
were operational. The lower number gives the total non-playground time which was searched for 
inspiral triggers. The shaded region corresponds to the data used in the S2 MACHO search. 
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Figure 41 : The plot shows the size of the template bank, generated with a minimal match of 97%, for 
various values of TOmin- 



-8/3 



As described in section 5.3.1 the size of the template bank is proportional 
to m^^" , where m^i^ is the mass parameter of the smallest equal mass binary in the template 
bank. Using these data, it was decided that the lowest mass accessible was mj^m — 0-2 Mq, given 
the available computational resources. 
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Figure 42 : The bottom right plot shows the variation in the size of the MACHO search template 
bank over the course of the S2 run. As described in the text, the template bank is generated using 
LI data to cover a region of parameter space from 0.2 to 1.0 Mq (component mass) at 95% 
minimal match. This template bank is used to filter the LI data in the triggered search pipeline. 
The other three plots show the variation in distance to which the three LIGO interferometers can 
see an optimally oriented binary at signal-to-noise ratio 8 over the S2 run. Since this is a function 
of the masses of the binary, this range is shown for a (0.1, 0.1) Mq and a (0.5, 0.5) Mq binary. 
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Figure 43 : The figure shows a histogram of all the triggers generated from the HI data using the 
triggered search (i.e. no coincidence with LI or H2 has been applied to the triggers). The signal-to- 
noise threshold is = 7 and the parameters of the veto are p — 15, S'^ = 0.04, S = 12.5, as in 
the S2 binary neutron star search. If the interferometer data is Gaussian, then wc would expect the 
histogram to be monotonically decreasing with increasing signal-to-noise ratio; however, there is a 
pronounced "hump" in the histogram at p « 9 suggesting some non-Gaussian feature in the data. 
By lowering the value of S to 5, we can remove this feature from the histogram, but we must be 
careful in doing so that we do not reduce the detection efficiency of the pipeline. 
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Figure 44 : The figure shows the results of a small Monte Carlo simulation used to test the detection 
efficiency of the pipeline using HI triggers (i.e. no coincidence with LI or H2 has been applied 
to the triggers). The signal-to-noise threshold is = 7 and the parameters of the veto are 
p — 15, (5^ — 0.04, S — 12.5. Found injections arc shown with a +, missed injections are shown 
with a X and the masses of the injection are shown as the x and y coordinates. We would expect 
to miss any injections with a mass component below 0.2 due to the coverage of the template 
bank; however injections in the region inside the bank should be detected, unless they are at an 
effective distance larger than the range of the interferometer. The missed injections that we would 
expect to find are color coded according to the effective distance at which they are injected. Several 
injections are missed as they are at a large effective distance (e.g. the injection at (0.48,0.42) Mq); 
however there are may missed injections at distances < 200 kpc which should be detectable in the 
HI data (e.g. the injection at (0.88,0.91) M©). Investigation of the missed injections showed they 
had large values of signal-to-noise ratio, but were vetoed by the test. This suggests that the 
search parameters used must be re-tuned to increase the detection efficiency. 
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Figure 45 : The plots in this figure show the values of p and S — x^/(15 + 5^p^) for the inspiral 
triggers corresponding to injected signals found by the triggered search pipeline. The upper plot 
shows LI triggers and the lower plot shows HI triggers (which correspond to the found injections 
of figure 44 1. No coincidence has been applied to the HI triggers at this stage; however they are 
generated using template banks produced from LI triggers. The color of each trigger shows the 
effective distance at which it was injected. Both plots are generated with a signal-to-noise threshold 
of p* = 7, and the parameters of the veto werep = 15, 5'^ = 0.04 and Sli = 12.5, Sli = 5.0, values 
chosen based on the tuning of the S2 binary neutron star search. It can be seen that, at a given 
signal-to-noise ratio, the HI triggers typically have higher values of S than the LI triggers. This 
is due to the a larger mismatch between the injected signal and the templates in the HI triggered 
bank. 




Figure 46 : The plots in this figure should show the observed values of p and S = / {lb + 5^ p^) for 
the inspiral triggers corresponding to injected signals found by the triggered search pipeline using 
the final set of parameters chosen. The upper plot shows LI triggers and the lower plot shows HI 
triggers. No coincidence has been applied to the HI triggers at this stage; however they are generated 
using template banks produced from LI triggers. These plots should be compared to those shown in 
figure 45 By tuning the value of 5^ to 0.2, it can be seen that much lower values of S are obtained 
for the HI injections. This suggests that we could further reduce the threshold S*, although this 
was not done as no coincident triggers were found in the playground data and a looser value of 5 
allowed us to probe the region slightly outside the template bank parameter space. 
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Parameter Description value 

/hp High Pass Filter Frequency 100 Hz 

Ohp High Pass Filter Order 100 Hz 

/low Low Frequency Cutoff 100 Hz 

JTimin Template bank lower component mass 0.2 Mq 

TTimax Template bank Upper Component mass 1.0 M© 

M LI template bank minimal match 0.95 

LI signal-to-noise ratio threshold 7.0 

LI veto threshold 3.1 

Phi HI signal-to-noise ratio threshold 7.0 

HI veto threshold 5.0 

Ph2 H2 signal-to-noise ratio threshold 7.0 

S^2 H2 veto threshold 10.0 

p Number of bins in x^ veto 15 

5^ X^ veto mismatch parameter 0.2 

Sm Trigger mass coincidence parameter 0.0 

6tiin H1-H2 trigger time coincidence parameter 0.001 s 

5tLH Ll-Hl, L1-H2 trigger time coincidence parameter 0.011 s 

khh H1-H2 trigger amplitude coincidence parameter 0.5 

Klh Ll-Hl, L1-H2 trigger amplitude coincidence parameter 1000.0 

e Trigger amplitude coincidence parameter 2.0 



Table 7 : A complete list of the parameters that were selected at the various stages of the pipeline. 
These values are justified in the text. 




Figure 47 : This figure shows the results of the Monte Carlo simulation used to measure the efhciency 
of the pipeline once parameter tuning had been completed; these detected triggers have survived 
all threshold and coincidence tests. The injections that are detected are shown as circles on the 
upper plot. The lower plot shows the injections that were not detected: stars correspond to missed 
injections in the triple coincident data, upward pointing triangles to the Ll-Hl data and downward 
pointing triangles to the L1-H2 data. The x and y coordinates are the mass parameters mi and TO2 
of each injection, respectively. The color of each injection represents the effective distance in the 
Hanford interferometers at which the waveform was injected (since the LHO interferometers limit 
the sensitivity of the search). The horizontal and vertical red lines show the edge of the template 
bank parameter space. 
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Figure 48 : The upper two plots show the efficiency of the pipeline e and the loss of the search 1 — e 
as a function of the injected signal chirp mass A4 measured by the Monte Carlo Simulation. The 
lower plot shows the chirp mass and effective distance in LHO of the injections used to measure the 
pipeline efficiency; detected injections are shown with a + and missed injections are shown with a x . 
It can be seen that the efficiency of the pipeline is unity very close to unity for high values of Ai and 
falls as the chirp mass decreases. There appears to be an anomalously large value of e at w 0.18, 
however it can be seen from the lower plot that there were comparatively few injections at this chirp 
mass, so this may be an effect of small number statistics. Further Monte Carlo simulations will be 
able to test this hypothesis. 
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Figure 49 : The panels in this figure compare the measured values of eff'ective distance, chirp mass 
and end time with the known values for the injected waveforms in the Monte Carlo simulation in 
the LI detector. 
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Figure 50 : The panels in this figure compare the measured values of effective distance, chirp mass 
and end time with the known values for the injected waveforms in the Monte Carlo simulation in 
the HI detector. 
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Figure 51 : The panels in this figure compare the measured values of effective distance, chirp mass 
and end time with the known values for the injected waveforms in the Monte Carlo simulation in 
the H2 detector. 
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Figure 52 : The signal-to-noise ratio of the background triggers produced by 20 time-sHdes. No triple 
coincident background triggers were observed. The colors are color coded depending on whether they 
were found in the triple, Ll-Hl double of Ll-Hl double triggers coincident data. No background 
triggers were found coincident in all three detectors, so the triggers from the triple coincident data 
set are from Ll-Hl coincidence only. 
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Figure 53 : The plots in this figure compare the signal-to-noise ratios of the background triggers to 
those of the triggers corresponding to software injections from the Monte Carlo simulation. Notice 
in the upper plot that the signal-to-noise ratio of detected injections in LI is a factor of ^ 2 higher 
than the signal-to-noise ratio in the LHO detectors, due to the greater sensitivity of LI. The lower 
plot shows a magnification of the low signal-to-noise ratio region. It can be seen that background 
triggers generally have a larger signal-to-noise ratio in the LHO detectors, suggesting the coherent 
statistic described in the text that gives greater weight to the LI signal-to- noise ratio. 
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Chapter 8 



Conclusion 



Although the upper hmit that we have placed on the rate of binary black hole MACHO 
inspirals in the galaxy is lower than the upper bound of the predicted rates, the LIGO 
interferometers were not at design sensitivity when the S2 data was taken. At present, 
the sensitivities of the instruments are significantly better than during S2, as can be 



seen from figure 54, and progress on reducing noise in the interferometers continues 
apace. The increase in detector sensitivity makes a larger volume of the Universe 
accessible to searches for binary inspirals. In addition to this, the amount of data is 
also increasing as the interferometers become more stable. 

These improvements in the instruments will increase the chance of detecting grav- 
itational waves from binary inspirals. If the rates of binary black hole MACHO 
coalescence are truly as high as predicted, then initial LIGO would stand an excel- 
lent chance of detecting an inspiral. The first detection of gravitational waves will 
be a major scientific breakthrough and will yield and enormous amount of scientific 
information, particularly if the detection came from a binary black hole MACHO. 
The length of binary black hole MACHO inspirals in the sensitive band of the in- 
terferometer will allow extremely accurate parameter estimation as well as tests of 
post-Newtonian theory. For systems with total mass greater than ~ 0.64 Mq LIGO 
will be sensitive to the coalescence of the binary and will be able to study the strong 
gravitational field effects when two binary black holes merge. When this is coupled 
with the accurate parameter estimation available from the earlier part of the wave- 
form, the inspiral of a binary black hole MACHO could be an excellent laboratory for 
General Relativity. A detection would also impact the studies of halo dark matter and 
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early universe physics, providing a MACHO component to the halo and suggesting 
that primordial black holes do indeed form in the universe. 

In the absence of detection, the improvements in detector sensitivity will dra- 
matically improve the upper limits placed on the rate of binary black hole MACHO 
inspirals. Once these rates are below the predicted rates, we may begin to use obser- 
vations from gravitational wave interferometers to constrain the fraction of galactic 
halos in the form of primordial black hole MACHOs. While this may not be as 
significant as a detection, it will still be of interest to the astrophysical community. 
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Figure 54 : Comparison of the best sensitivities of the LIGO interferometers between science runs. 
The sohd curve shows the design sensitivity for the 4 km interferometers: the LHO 4 km is only a 
factor of ~ 2 away from design at 100 Hz during S3. 
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